<repo_name>AlexChuk/1D_HF_Discharge<file_sep>/Read_mesh.h
//Защита от повторного включения заголовка
#ifndef READ_MESH_H
#define READ_MESH_H

# define Gf_P 10  //geom_factors(Poisson)

void read_mesh(int);
void calc_geom_rz(int);
void calc_geom_xy(int)

#endif
<file_sep>/Gas_calc.cpp
# include "1D_MainFun.h"
# include "Gas_calc.h"

double HCpSi[3][Nmax];

void gas_HCpSi_calc(double Tgas,int N)//нахождение внутренней энергии,энтальпии и теплоёмкостей по зад. температуре из аппроксиммаций Cp(T)
{
	double Hi,Cpi,Si,Mm,Cpt;
	int l,n,x;

	//T0 = Temp0;
	for(n=1;n<N;n++)
	{
		Mm = Mi[n];//[г]
		//Расчёт характеристик в точке T0
		if(Tgas>1000)
			x=0;
		else
			x=1;

		Cpi = 0;
		Hi = 0;
		Si = 0;
		for(l=0;l<5;l++)
		{
			Cpt = CXi[n][x][l]*pow(Tgas,l);

			Cpi += Cpt;//безразм
			Hi += Cpt/(l+1);//безразм

			if(l!=0)
				Si += Cpt/l;//безразм
		}
		Cpi = Cpi*kb/Mm;//[эрг/г*К]
		Hi = (Hi*kb*Tgas+CXi[n][x][5]*kb+CXi[n][x][7]*1.6e-12)/Mm;//[эрг/г]
		Si = (Si+CXi[n][x][0]*log(Tgas)+CXi[n][x][6])*kb/Mm;//[эрг/г*К]

		HCpSi[0][n] = Hi;//[эрг/г]
		HCpSi[1][n] = Cpi;//[эрг/г*К]
		HCpSi[2][n] = Si;//[эрг/г*К]
	}
}
void gas_TP_calc(double *Nn,double *Xn,int N,double *Pgas,double *Tgas,double *dTgas,double *Ngas,double *Rogas,double *Hgas,double *Nel,double *dNel,double *Tv)//расчёт температуры газа((H,P)-const)
{
	int n;
	double Xi[N],Nin,Roin,Tin,Pin,Hin,XMi,Xmi[N],hin;
	double Nout,Pout,Tout,Roout,Hout;
    double Hi,Cpi;
	double ftn,Ftn,Tnn,Tn,dT;

	*dNel = *Nel;

	//Присвоение значений из адресов переменных
	Tin = *Tgas;
	Pin = *Pgas;
	Hin = *Hgas;

	Nin = 0.0;//*Ngas;//0.0;
	for(n=0;n<N;n++)
		Nin += *(Nn+n*(LEN+2));//Nin += Ni[n*(LEN+2)];

	XMi = 0.0;
    for(n=0;n<N;n++)
    {
        Xi[n] = *(Nn+n*(LEN+2))/Nin;//*(Xn+n*(LEN+2));////*(Xn+n*(LEN+2));//*(Nn+n*(LEN+2))/Nin;//*(Xn+n*(LEN+2));//Xi[n] = Ni[n*(LEN+2)]/Nin;
        Xmi[n] = Xi[n]*Mi[n];
        XMi += Xmi[n];
    }
    Roin = XMi*Nin;

    //Расчёт температуры_методом Ньютона
	Tn = Tin;//300;
	//hin = Hin/Nin;
	Tnn = Tin;
	do
	{
		if(!Tnn==0)
            Tn = Tnn;

		Ftn = 0;
		ftn = 0;
		gas_HCpSi_calc(Tn,N);
		for(n=1;n<N;n++)//Npos+Nneg
		{
			Hi = HCpSi[0][n];//[эрг/г]
			Cpi = HCpSi[1][n];//[эрг/г*К]

			Ftn += Xmi[n]*Nin*Hi;//Roi1[n]*Hi;//Xi[n]*Hi;//
			ftn += Xmi[n]*Nin*Cpi;//Roi1[n]*Cpi;//Xi[n]*Cpi;//
		}

		//Ftn = Ftn*Nn - Hin;
		Ftn = Ftn - Hin;//hin;//

		//Tnn = Tn - Ftn/(ftn*Nn);
		Tnn = Tn - Ftn/ftn;

		dT = fabs(Tnn-Tn);

	}while(dT>0.001);//0.001
	Tout = Tnn;//Tin//300;
	if(Tout<300.0)
        Tout = 300.0;

	//Isobaric process**************************************
	Pout = Pin;

	//Return to using variables:
	Nout = Pout/(kb*Tout);
	//Nout = Nin;
    Roout = XMi*Nout;

	gas_HCpSi_calc(Tout,N);
	Hout = 0;
	for(n=1;n<N;n++)
		Hout += Xmi[n]*HCpSi[0][n];//[эрг/cm^3]
    Hout *= Nout;
    //Hout = Hin;//*Nout/Nin;

	//Return to using variables:
	*Pgas = Pout;
	*Tgas = Tout;
	*Ngas = Nout;
	*Rogas = Roout;
	*Hgas = Hout;
    *dTgas = fabs(Tin-Tout);

	for(n=0;n<N;n++)
    {
        *(Nn+n*(LEN+2)) = Nout*Xi[n];//Ni[n*(LEN+2)] = Nout*Xi[n];
        *(Xn+n*(LEN+2)) = Xi[n];
    }

    //концентрация электронов
    *Nel = Nout*Xi[0];
    *dNel = fabs(*dNel-*Nel);

    //Tv-calculation
    /*for(n=0;n<N;n++)
    {
        if(!strcmp(Spec[n],"N2(0)"))
            break;
    }*/
    n=11;
    *Tv = fabs(HCpSi[0][n+1]-HCpSi[0][n])*Mi[n]/kb/log((*(Nn+n*(LEN+2)))/(*(Nn+(n+1)*(LEN+2))));//*Tv = fabs(HCpSi[0][n+1]-HCpSi[0][n])*Mi[n]/kb/log(Ni[n*(LEN+2)]/Ni[(n+1)*(LEN+2)]);

	//************************************************************
}
void gas_LenPrint(double *Xi,double *Ni,int N,int Npos,int Nneg,double *Pgas,double *Tgas,double *Ngas,double *Hgas,double *Rogas,double *Nel,double *Te,double *Tv,double *Er,double *Fir,double *Jel,double tic)//запись в файл
{
	int i,k,n,nmax;
	FILE *log;
	int I[] = {1,int(LEN/2),LEN};
	double Ch;
	char ChL[][5] = {"_0","_L/2","_L"};

	nmax = N-1;//12;

	///Veusz_data_printing*****************************************************
	log = fopen("GUIData.txt", "w");

    fprintf(log,"r,cm\t");
    for(n=0;n<=nmax;n++)
        fprintf(log,"%s\t",Spec[n]);
    fprintf(log,"Charge\n");

    ///Concentrations
    for(i=0;i<LEN+2;i++)
    {
        //fprintf(log,"%.3lf\t",(i+0.5)*Len/(LEN+1));//old
        fprintf(log,"%.3lf\t",0.5*(l[i+1]+l[i]));
        for(n=0;n<=nmax;n++)
            fprintf(log,"%.2e\t",Ni[n*(LEN+2)+i]);

        Ch = -Ni[i] ;//electrons
        for(n=1;n<=Npos;n++)//positive ions
            Ch += Ni[n*(LEN+2)+i];
        for(n>Npos;n<=Npos+Nneg;n++)//negative ions
            Ch += -Ni[n*(LEN+2)+i] ;

        //fprintf(log,"%.2e\n",Ch);
        fprintf(log,"%.2e\n",fabs(Ch));
    }
    fprintf(log,"\n");

    ///Mole_fractions
    for(n=0;n<=nmax;n++)
        fprintf(log,"x(%s)\t",Spec[n]);
    fprintf(log,"\n");

    for(i=0;i<LEN+2;i++)
    {
        for(n=0;n<=nmax;n++)
            fprintf(log,"%.2e\t",Xi[n*(LEN+2)+i]);
        fprintf(log,"\n");
    }
    fprintf(log,"\n");

    //fprintf(log,"Tgas,K\tTv,K\tTe,K\tTe,eV\tE/N,Td\n");
    fprintf(log,"Ngas,cm-3\tTgas,K\tTe,K\tTe,eV\n");
    for(i=0;i<LEN+2;i++)
        //fprintf(log,"%.2e\t%.2lf\t%.1lf\t%.1lf\t%.3lf\t%.2lf\n",Ngas[i],Tgas[i],Tv[i],Te[i]*eV_K,Te[i],E*1.0e17*Eabs/Ngas[i]);
        fprintf(log,"%.2e\t%.2lf\t%.1lf\t%.3lf\n",Ngas[i],Tgas[i],Te[i]*eV_K,Te[i]);
    fprintf(log,"\n");

    fprintf(log,"Fi,V\tEr/N,Td\tEr,V/cm\t");
    for(i=0;i<LEN+1;i++)
        fprintf(log,"%.3e\t%.3e\t%.3e\n",Fir[i]*Eabs,Er[i]*Eabs*1.e17/Ngas[i],Er[i]*Eabs);
    fprintf(log,"%.3e\n",Fir[LEN+1]*Eabs);
    fprintf(log,"\n");

    /*//VDF_print
    fprintf(log,"n\t");
    for(i=0;i<3;i++)
        //fprintf(log,"VDF(%.2lfcm)\t",l[I[i]]);
        fprintf(log,"VDF%s\t",ChL[i],l[I[i]]);
    fprintf(log,"\n");
    for(n=11;n<N;n++)
    {
        fprintf(log,"%d\t",n-11);
        for(i=0;i<3;i++)
            fprintf(log,"%.2e\t",Ni[n*(LEN+2)+I[i]]/Ni[11*(LEN+2)+I[i]]);
        fprintf(log,"\n");
    }
    fprintf(log,"\n");*/

    fclose(log);
}
void gas_TimePrint(double *Ni,int N,double Pgas,double Tgas,double Ngas,double Te,double Tv,double Jel,double Icalc,double Vext,double Wext,double Pcalc,double PowTav,double IcalcTav,double MNT,double MNTr,double tic,double etic,double dte,double dt)//запись в файл
{
    int i,k,n;
    double NRF = 1.0;
	FILE *log;

    ///Output_on_monitor*******************************************************
	printf("Time = %.2e[s]\n",tic);
	printf("LENaveraged_Data:\n");
	printf("P = %.1lf[Torr]\tT = %.1lf[K]\tTv = %.1lf[K]\n",Pgas/p0,Tgas,Tv);
	printf("I = %.2lf[mA]\tVext = %.2lf[kV]\tPower = %.2e[W/cm^2]\n",Icalc*1e3,Vext*Eabs/1.e3,Pcalc);
	printf("Xe = %.2e\t%s = %.2e[cm-3]\t%s = %.2e[cm-3]\t%s = %.2e[cm-3]\n\n",Ni[0]/Ngas,Spec[0],Ni[0],Spec[1],Ni[1],Spec[4],Ni[4]);

    ///Output_to_file**********************************************************
	if(tic==0.0)
    {
        log = fopen("Gas_TimeData.txt", "w");

        //fprintf(log,"t,s\tdte,s\tdt,s\tPgas(t),Torr\tTgas(t),K\tTv(t),[K]\tTe(t),[K]\tEz(t),V/cm\tE/N(t),Td\tJel(t),[A/cm2]\tIel(t),[mA]\tNgas(t),[cm-3]\tXel\tMNT,[g*cm-3*K]\tMNTr,[g*cm-3*K]\t");
        fprintf(log,"t,s\tt_el,s\tdte,s\tdt,s\tNRF\tPgas(t),Torr\tTgas(t),K\tTe(t),[K]\tJel(t),[A/cm2]\tIel(t),[mA]\tVext(t),[kV]\tPower(t),[W/cm^2]\t<Iel>RF,[mA]\t<Power>RF,[W/cm^2]\tNgas(t),[cm-3]\tXel\tMNT,[g*cm-3*K]\tMNTr,[g*cm-3*K]\t");///w/o_Tv
        for(n=0;n<N;n++)//15
            fprintf(log,"%s(t)\t",Spec[n]);
        fprintf(log,"\n");
    }
    else
        log = fopen("Gas_TimeData.txt", "a+");

    if(Wext!=0)
    {
        Vext *= sin(2*pi*Wext*etic);
        NRF = etic*Wext;
    }

    //fprintf(log,"%.7e\t%.2e\t%.2e\t%.1lf\t%.1lf\t%.1lf\t%.1lf\t%.1lf\t%.1lf\t%.2e\t%.1lf\t%.2e\t%.2e\t%.5e\t%.5e\t",tic,dte,dt,Pgas/p0,Tgas,Tv,Te*eV_K,E*Eabs,E*1.0e17*Eabs/Ngas,Jel*eKl/e,Icalc*1e3,Ngas,Nel/Ngas,MNT,MNTr);
    fprintf(log,"%.7e\t%.7e\t%.2e\t%.2e\t%.2lf\t%.1lf\t%.1lf\t%.1lf\t%.2e\t%.3lf\t%.3lf\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.5e\t%.5e\t",tic,etic,dte,dt,NRF,Pgas/p0,Tgas,Te*eV_K,Jel*eKl/e,Icalc*1.e3,Vext*Eabs/1.e3,Pcalc,IcalcTav*1.e3,PowTav,Ngas,Ni[0]/Ngas,MNT,MNTr);///w/o_Tv
    for(n=0;n<N;n++)//15
        fprintf(log,"%.2e\t",Ni[n]);
    fprintf(log,"\n");

	fclose(log);
}
void gas_SavePrint(double *Ni,int N,double *Pgas,double *Tgas,double *Nel,double *Te,double *Fir,double tic)//запись в save-файл
{
	int i,k,n;
	FILE *save;

    //запись параметров газа*******************************************************************
	save = fopen("Save_data.txt", "w");

	fprintf(save,"Time,sec\t%.2e\n",tic);

    //Data_print
    fprintf(save,"Pgas,Torr\t%.1lf\n",Pgas[i]/p0);

    fprintf(save,"Tgas,K\t");
    for(i=0;i<LEN+2;i++)
        fprintf(save,"%.1lf\t",Tgas[i]);
    fprintf(save,"\n");

    fprintf(save,"Te,eV\t");
    for(i=0;i<LEN+2;i++)
        fprintf(save,"%.2lf\t",Te[i]);
    fprintf(save,"\n");

    //Field
    fprintf(save,"Fi,V\t");
    for(i=0;i<LEN+2;i++)
        fprintf(save,"%.2e\t",Fir[i]);//(LEN+2) - due to 2 additional virtual points
    fprintf(save,"\n");

    //Ni[n]
    for(n=0;n<N;n++)
    {
        fprintf(save,"%s\t",Spec[n]);//,cm-3
        for(i=0;i<LEN+2;i++)
            fprintf(save,"%.2e\t",Ni[n*(LEN+2)+i]);//(LEN+2) - due to 2 additional virtual points
        fprintf(save,"\n");
    }

    //electron_data

    /*//VDF_print
    int I[] = {1,int(LEN/2),LEN};
    for(i=0;i<3;i++)
    {
        fprintf(save,"VDF(l=%.2lfcm)\t",l[I[i]]);
        for(n=11;n<N;n++)
            fprintf(save,"%.2e\t",Ni[n*(LEN+2)+I[i]]/Ni[11*(LEN+2)+I[i]]);
        fprintf(save,"\n");
    }
    fprintf(save,"\n");

    fprintf(save,"\n\n");*/

	fclose(save);
}
void gas_LenAverage(int N,char *Geom,double *Xi,double *Ni,double *Ngas,double *Tgas,double *Te,double *Tv,double *Xi_av,double *Ni_av,double *Ng_av,double *Tg_av,double *Te_av,double *Tv_av)
{
    int n;
    for(n=0;n<N;n++)
    {
        Xi_av[n] = LENaverage(&Xi[n*(LEN+2)],Geom);
        Ni_av[n] = LENaverage(&Ni[n*(LEN+2)],Geom);
    }

    *Ng_av = LENaverage(Ngas,Geom);
    *Tg_av = LENaverage(Tgas,Geom);
    *Te_av = LENaverage(Te,Geom);
    *Tv_av = LENaverage(Tv,Geom);
}


<file_sep>/Chemistry_calc.cpp
# include "1D_MainFun.h"
# include "Gas_calc.h"
# include "Chemistry_calc.h"
# include "1Dmesh.h"

int VMax[Nmax];
int ChNRe[NRmax][3];//массив номеров реагентов хим.реакций
int ChM[Nmax][NRmax],Chem[Nmax][NRmax+3];//массив мультипликаторов реакции,массив необходимых реакций
int ChRad[NRmax];//массив номеров радиационных реакций
double dEch[NRmax];//массив дефектов энергий реакций
double te_ch = 1.0, tg_ch = 1.0, tic_ch;
int Iech,Igch;

char Rtype[NRmax][10];//массив типов хим.реакций

double Kchi[NRmax][8];//массив скоростей,констант и коэффициентов скоростей хим.реакций
double RP[Nmax][LEN+2][NRmax],RC[Nmax][LEN+2][NRmax];//,RPp[LEN][Nmax][NRmax],RCp[LEN][Nmax][NRmax];
char RName[NRmax][100];

int chem_make_react(int Nedf,char *ChemFile)//формирование общего файла хим. процессов
{
	int j,n,nr,nRt,Rt,rev,J,Nn,Nvib,Nchem,k,Ki;
	char CHstr[20],Cmt[100],RTypCmt[10][100],Tcmt[10];
	char symb[] = "------------------------------------------------------------";
	char Creac[10][10],Micmt[10][10];
	double x1[10],CKj[10];

	FILE *react;
	//react = fopen("N2_chemistry.txt", "r");
	//react = fopen("N2_chemistry_new.txt", "r");
	//react = fopen("N2_chemistry_new-2.txt", "r");
	//react = fopen("N2_chemistry_new-test.txt", "r");
	//react = fopen("N2_chemistry_new-test.txt", "r");
	//react = fopen("N2_chemistry_old.txt", "r");
	//react = fopen("test-react.txt", "r");
	react = fopen(ChemFile, "r");

	FILE *rset;
	rset = fopen("ReactionSet.txt", "a+");

	//Printing diagnosic message************************************************
	printf("Chemistry info:\n\n");

	//считывание комментария
	fscanf(react,"%s",&Cmt);
	do
	{
		fscanf(react,"%s",&Cmt);
	}while(strcmp(Cmt,symb)!=0);

	Rt = 0;
	Nchem = Nedf;
	Nvib = 0;
	fscanf(react,"%s",&Cmt);
	while(strcmp(Cmt,"END")!=0)
	{
		strcpy(RTypCmt[Rt],Cmt);//RTypCmt[Rt] = Cmt;

		//printing name of reactions' block
		fprintf(rset,"\n%s\n",symb);
		fprintf(rset,"\n%s\n",RTypCmt[Rt]);
		fprintf(rset,"\n%s\n",symb);

		fscanf(react,"%s",&Cmt);
		nRt = 0;
		fscanf(react,"%s",&CHstr);
		do//read reactions inside the block
		{
			//read i-th reaction
			//left part of the reaction
			j = 0;
			while((strcmp(CHstr,"->")!=0)&&(strcmp(CHstr,"<->")!=0))
			{
				strcpy(Creac[j],CHstr);//Creac[j] = CHstr;
				j += 1;
				fscanf(react,"%s",&CHstr);
			}

			//reverse or not?
			if(!strcmp(CHstr,"<->"))
				rev = 1;
			else
				rev = 0;

			strcpy(Creac[j],"->");//Creac[j] = "->";
			j += 1;

			//right part of the reaction
			fscanf(react,"%s",&CHstr);
			while(strcmp(CHstr,";")!=0)
			{
				strcpy(Creac[j],CHstr);//Creac[j] = CHstr;
				j += 1;
				fscanf(react,"%s",&CHstr);
			}
			J = j;

			//coefficient
			fscanf(react,"%s",&Tcmt);
            if(!strcmp(Tcmt,"Te"))
                Ki = 3;
            if(!strcmp(Tcmt,"Tgas") || !strcmp(Tcmt,"Radiat"))
                Ki = 3;
            if(!strcmp(Tcmt,"VV"))
                Ki = 3;
            if(!strcmp(Tcmt,"VV'"))
                Ki = 3;
            if(!strcmp(Tcmt,"VT"))
                Ki = 3;

            for(k=0;k<Ki;k++)
                fscanf(react,"%lf",&CKj[k]);

			//different components
			fscanf(react,"%s",&CHstr);
			n = 0;
			Nn = 0;
			if(!strcmp(CHstr,"M:"))
			{
				fscanf(react,"%s",&CHstr);
				while(strcmp(CHstr,";")!=0)
				{
					strcpy(Micmt[n],CHstr);//Micmt[n] = CHstr;
					fscanf(react,"%lf",&x1[n]);
					n += 1;

					fscanf(react,"%s",&CHstr);
				}

				fscanf(react,"%s",&CHstr);
				Nn = n-1;
			}

			//add reactions to the ReactionSet.txt
            nr = Nchem+nRt;
            if((!strcmp(Tcmt,"VV"))||(!strcmp(Tcmt,"VV'"))||(!strcmp(Tcmt,"VT")))
            {
                fclose(rset);

                Nvib = chem_VV_VT_make(nr,Tcmt,Creac[0],Creac[2],CKj,Ki);

                nRt += Nvib;
            }
            else
            {
                for(n=0;n<=Nn;n++)
                {
                    nr += 1;
                    fprintf(rset,"R%d\t",nr);
                    for(j=0;j<J;j++)
                    {
                        if((!strcmp(Creac[j],"M"))&&(Nn!=0))
                            fprintf(rset,"%s\t",Micmt[n]);
                        else
                            fprintf(rset,"%s\t",Creac[j]);
                    }
                    fprintf(rset,";\t");

                    if(Nn==0)
                        fprintf(rset,"%s\t%.3e\t",Tcmt,CKj[0]);
                    else
                        fprintf(rset,"%s\t%.3e\t",Tcmt,CKj[0]*x1[n]);
                    for(k=1;k<Ki;k++)
                        fprintf(rset,"%.3lf\t",CKj[k]);
                    fprintf(rset,"%s\n",CHstr);
                    nRt += 1;

                    if(rev==1)
                    {
                        nr += 1;
                        fprintf(rset,"R%d\t",nr);
                        for(j=J-1;j>=0;j--)
                        {
                            if((!strcmp(Creac[j],"M"))&&(Nn!=0))
                                fprintf(rset,"%s\t",Micmt[n]);
                            else
                            fprintf(rset,"%s\t",Creac[j]);
                        }
                        fprintf(rset,";\t");
                        fprintf(rset,"Rev\t\t%s\n",CHstr);
                        nRt += 1;
                    }
                }
			}

			//!!***************************
			//введение исключения чтения
			//!!***************************

			fscanf(react,"%s",&Cmt);
			if(strcmp(Cmt,symb))
				strcpy(CHstr,Cmt);

		}while(strcmp(Cmt,symb)!=0);

		printf("%d\t%s were added to ReactionSet.txt\n",nRt,RTypCmt[Rt]);
		Rt += 1;
		Nchem += nRt;

		fscanf(react,"%s",&Cmt);
	}
	printf("Total of %d reactions were added to ReactionSet.txt\n\n%s\n\n",Nchem,symb);

	fclose(react);
	fclose(rset);

	return Nchem;
}
int chem_VV_VT_make(int Nchem,char Cmt[],char VSpec1[],char VSpec2[],double CKi[],int Ki)//формирование VV_VT_реакций
{
    int Nvib = 0,v1,v2,Vmax1,Vmax2,n,n1,n2,i,quant;
    char str[10],str1[10];

    FILE *rset;
	rset = fopen("ReactionSet.txt", "a+");

    //определение 1-го компонента и максимума колебательного номера
    for(n=0;n<Nmax;n++)
    {
        if(!strcmp(VSpec1,Spec[n]))
        {
            n1 = n;
            if(!strcmp(Cmt,"VT"))
                n1 = n-1;

            i = 0;
            memset(str, 0, sizeof(str));
            memset(str1, 0, sizeof(str1));
            while(VSpec1[i]!='(')
            {
                str[i] = VSpec1[i];
                i++;
            }

            Vmax1 = 0;
            do
            {
                Vmax1 += 1;
                sprintf(str1,"%s(%d)",str,Vmax1);
            }while(!strcmp(str1,Spec[n1+Vmax1]));

            VMax[n1] = Vmax1;

            break;
        }
    }

    //определение 2-го компонента
    for(n=0;n<Nmax;n++)
    {
        if(!strcmp(VSpec2,Spec[n]))
        {
            n2 = n;
            break;
        }
    }

	if(!strcmp(Cmt,"VV"))
    {
        quant = n2 - n1;

        for(v1=0;v1<Vmax1-quant;v1++)//(v1=0;v1<Vmax1-1;v1++)
        {
            for(v2=quant;v2<Vmax1;v2++)//(v2=1;v2<Vmax1;v2++)
            {
                fprintf(rset,"R%d\t %s\t +\t %s\t->\t %s\t +\t %s\t;\t",Nchem+1+Nvib,Spec[n1+v1],Spec[n1+v2],Spec[n1+v1+1],Spec[n1+v2-1]);
                if(Nvib==0)
                {
                    fprintf(rset,"VV\t %.2e\t",CKi[0]);
                    for(i=1;i<Ki;i++)
                        fprintf(rset,"%.2lf\t",CKi[i]);
                }
                else
                    fprintf(rset,"None\t");
                fprintf(rset,"//ref\n");

                Nvib += 1;
            }
        }
    }

    if(!strcmp(Cmt,"VV'"))
    {
        i = 0;
        memset(str, 0, sizeof(str));
        memset(str1, 0, sizeof(str1));
        while(VSpec2[i]!='(')
        {
            str[i] = VSpec2[i];
            i++;
        }

        n2 -= 1;
        Vmax2 = 0;
        do
        {
            Vmax2 += 1;
            sprintf(str1,"%s(%d)",str,Vmax2);
        }while(!strcmp(str1,Spec[n2+Vmax2]));

        VMax[n2] = Vmax2;

        for(v1=0;v1<Vmax1-1;v1++)
        {
            for(v2=1;v2<Vmax2;v2++)
            {
                fprintf(rset,"R%d\t %s\t +\t %s\t->\t %s\t +\t %s\t;\t",Nchem+1+Nvib,Spec[n1+v1],Spec[n2+v2],Spec[n1+v1+1],Spec[n2+v2-1]);
                if(Nvib==0)
                {
                    fprintf(rset,"VV'\t %.2e\t",CKi[0]);
                    for(i=1;i<Ki;i++)
                        fprintf(rset,"%.2lf\t",CKi[i]);
                }
                else
                    fprintf(rset,"None\t");
                fprintf(rset,"//ref\n");

                Nvib += 1;

                fprintf(rset,"R%d\t %s\t +\t %s\t->\t %s\t +\t %s\t;\tNone\t",Nchem+1+Nvib,Spec[n1+v1+1],Spec[n2+v2-1],Spec[n1+v1],Spec[n2+v2]);
                fprintf(rset,"//ref\n");

                Nvib += 1;
            }
        }
    }

    if(!strcmp(Cmt,"VT"))
    {
        for(v1=1;v1<Vmax1;v1++)
        {
            fprintf(rset,"R%d\t %s\t +\t %s\t->\t %s\t +\t %s\t;\t",Nchem+1+Nvib,Spec[n1+v1],Spec[n2],Spec[n1+v1-1],Spec[n2]);
            if(Nvib==0)
            {
                fprintf(rset,"VT\t %.2e\t",CKi[0]);
                for(i=1;i<Ki;i++)
                    fprintf(rset,"%.2lf\t",CKi[i]);
            }
            else
                fprintf(rset,"None\t");
            fprintf(rset,"//ref\n");

            Nvib += 1;

            fprintf(rset,"R%d\t %s\t +\t %s\t->\t %s\t +\t %s\t;\tNone\t",Nchem+1+Nvib,Spec[n1+v1-1],Spec[n2],Spec[n1+v1],Spec[n2]);
            fprintf(rset,"//ref\n");

            Nvib += 1;
        }
    }

    fclose(rset);

    return Nvib;
}
void chem_read_react(int Nchem,int N)//считывание процессов
{
	int i,II,j,k,n,err,S_err,Num,n_rad;
	char Cmt[300],Rstr[10];
	char symb[] = "------------------------------------------------------------";

	FILE *react;
	react = fopen("ReactionSet.txt", "r");

	//считывание комментария
	fscanf(react,"%s",&Cmt);
	do
	{
		fscanf(react,"%s",&Cmt);
	}while(strcmp(Cmt,symb)!=0);

	//Reading of the reactions**********************************************************************
	Num = 0;
	S_err = 0;
	n_rad = 0;
	for(i=0;i<Nchem;i++)
	{
		if(i==0)
			fscanf(react,"%s",&Cmt);//номер реакции

		memset(RName[i], 0, sizeof(RName[i]));
		sprintf(RName[i],"(R%d)",i+1);

		//считывание левой части реакции+формирование матрицы левых частей
		j = 0;
        fscanf(react,"%s",&Rstr);
		while(strcmp(Rstr,"->")!=0)
		{
			strcat(RName[i],Rstr);

			if(!strcmp(Rstr,"M"))
			{
				ChNRe[i][j] = -1;
				j += 1;
			}
			else if(!strcmp(Rstr,"+"))
			{}
			else
			{
				err = 0;
				for(n=0;n<N;n++)
				{
					if(!strcmp(Rstr,Spec[n]))
					{
						ChNRe[i][j] = n+1;
						ChM[n][i] -= 1;
						err = 1;
						break;
					}
				}
				if(err==0)
				{
					S_err += 1;
					printf("Unexpected specie: Left - R#=%d\t %s\n",i+1,Rstr);
				}

				j += 1;
			}
			fscanf(react,"%s",&Rstr);
		}

        strcat(RName[i],"->");

		//считывание правой части реакции
		fscanf(react,"%s",&Rstr);
		while(strcmp(Rstr,";")!=0)
		{
			strcat(RName[i],Rstr);

			if(!strcmp(Rstr,"M"))
			{}
			else if(!strcmp(Rstr,"+"))
			{}
			else
			{
				err = 0;
				for(n=0;n<N;n++)
				{
					if(!strcmp(Rstr,Spec[n]))
					{
						ChM[n][i] += 1;
						err = 1;
						break;
					}
				}
				if(err==0)
				{
					S_err += 1;
					printf("Unexpected specie: Right - R#=%d\t %s\n",i+1,Rstr);
				}
			}

			fscanf(react,"%s",&Rstr);
		}

		fscanf(react,"%s",&Rtype[i]);//тип реакции
		if(!strcmp(Rtype[i],"EEDF"))
			II = 0;
		if(!strcmp(Rtype[i],"Te"))
			II = 3;
		if(!strcmp(Rtype[i],"Tgas"))
			II = 3;
        if(!strcmp(Rtype[i],"Radiat"))
        {
            II = 3;

            n_rad += 1;
            ChRad[0] = n_rad;
            ChRad[n_rad] = i;
        }
		if(!strcmp(Rtype[i],"Tgas-VT"))
			II = 3;
        if((!strcmp(Rtype[i],"VV"))||(!strcmp(Rtype[i],"VV'"))||(!strcmp(Rtype[i],"VT")))
            II = 3;
		if((!strcmp(Rtype[i],"Rev"))||(!strcmp(Rtype[i],"None")))
			II = 0;

		//считывание коэф-та скорости реакции
		for(k=0;k<II;k++)
			fscanf(react,"%lf",&Kchi[i][k]);
		fscanf(react,"%s",&Cmt);

		Num += 1;

		//считывание комментария
		fscanf(react,"%s",&Cmt);

		if(!strcmp(Cmt,symb))
		{
			do
			{
				fscanf(react,"%s",&Cmt);
			}while(strcmp(Cmt,symb)!=0);

			fscanf(react,"%s",&Cmt);
		}

	}
	fclose(react);

	//печать сообщения об ошибках или успехе
	if((S_err==0) && (Num==Nchem))
		printf("All reactions were carefully read\n\n");
	else
	{
		char a;

		printf("!!! Please, upgrade initial file!!Mistakes Sum = %d!!!\n",S_err);
		//printf("!!! Reactions read = %d, but written = %d!!!\n",Num,Nchem);

        scanf("%c",&a);
        exit(1);
	}

	printf("%s\n",symb);

	///Расчёт_дефектов_энергии_реакции*****************************
	gas_HCpSi_calc(300,N);
	for(i=0;i<Nchem;i++)
    {
        dEch[i] = 0.0;
        for(n=1;n<N;n++)
        {
            if(ChM[n][i]!=0)
                dEch[i] += -ChM[n][i]*HCpSi[0][n]*Mi[n];///[erg]//CXi[n][0][7]*1.6e-12;
        }
    }

	///Selection of necessary reactions****************************************************
	int I,Jr,Jc,nR,m;
	FILE *rrp;
	FILE *rcp;
	char str1[20];

	for(n=0;n<N;n++)
	{
		I = 0;
		Jr = 0;
		Jc = 0;
		for(i=0;i<Nchem;i++)
		{
			if(ChM[n][i]!=0)
			{
				Chem[n][I] = i;
				I += 1;

				if(ChM[n][i]>0)
					Jr += 1;
				else
					Jc += 1;
			}
		}
		Chem[n][Nchem] = I;//число реакций, в которых участвует компонент
		Chem[n][Nchem+1] = Jr;//число реакций(приход), в которых участвует компонент
		Chem[n][Nchem+2] = Jc;//число реакций(расход), в которых участвует компонент


		//Создание файлов скоростей интересующих компонент смеси**************************
		for(nR=0;nR<NR;nR++)
		{
			if(!strcmp(Spec[n],Spec_R[nR]))
			{
				//sprintf(str1, "%s%s", &Spec_R[nR], "_R+.txt");
				//rr = fopen(str1, "a+");
				sprintf(str1, "%s%s", Spec_R[nR], "_R(t)+.txt");
				rrp = fopen(str1, "w");
				//sprintf(str1, "%s%s", &Spec_R[nR], "_R-.txt");
				//rc = fopen(str1, "a+");
				sprintf(str1, "%s%s", Spec_R[nR], "_R(t)-.txt");
				rcp = fopen(str1, "w");

				//fprintf(rr,"t,s\t");
				fprintf(rrp,"tch,s\t");
				//fprintf(rc,"t,s\t");
				fprintf(rcp,"tch,s\t");

				for(j=0;j<I;j++)
				{
					m = Chem[n][j];

					if(ChM[n][m]<0)
					{
					//	fprintf(rc,"R%d\t",m+1);
						//fprintf(rcp,"R%d\t",m+1);
						fprintf(rcp,"%s\t",RName[m]);
					}
					else
					{
					//	fprintf(rr,"R%d\t",m+1);
						//fprintf(rrp,"R%d\t",m+1);
						fprintf(rrp,"%s\t",RName[m]);
					}
				}
				//fprintf(rr,"SUM\n");
				fprintf(rrp,"RSUM\n");
				//fprintf(rc,"SUM\n");
				fprintf(rcp,"RSUM\n");

				//fclose(rr);
				fclose(rrp);
				//fclose(rc);
				fclose(rcp);

				break;
			}
		}

	}

	//Logging*******************************************************************************
	FILE *log;
	log = fopen("Log_Chemistry.txt", "w");

	//запись компонент
	fprintf(log,"Components of the mixture:\n");
	for(n=0;n<N;n++)
		fprintf(log,"%d\t%s\n",n+1,Spec[n]);
	fprintf(log,"\n");

	//запись матрицы мультипликаторов реакции
	fprintf(log,"Matrix of multipliers:\n");
	for(n=0;n<N;n++)
	{
		fprintf(log,"%s\t\t",Spec[n]);
		for(i=0;i<Nchem;i++)
			fprintf(log,"%d\t",ChM[n][i]);
		fprintf(log,"\n");
	}
	fprintf(log,"\n");

	//запись матрицы необходимых реакций
	fprintf(log,"Matrix of necessary reactions:\n");
	for(n=0;n<N;n++)
	{
		fprintf(log,"%5s\t",Spec[n]);
		for(i=0;i<=Nchem+2;i++)
			fprintf(log,"%d\t",Chem[n][i]);
		fprintf(log,"\n");
	}
	fprintf(log,"\n");

	fclose(log);
}
void chem_const(double *Kch,int Nedf,int Nchem,int N,double Tel,double Tch,double tic)// расчёт констант скоростей реакций
{
	int j,n,Sum;
	double Kf;
	double Nv,kv10;

	int xRev=0;//счётчик запуска счета констант реакций с ID - Rev
    Tel *= eV_K;//[K]

	for(j=Nedf;j<Nchem;j++)
	{
		if(!strcmp(Rtype[j],"Te"))//запись констант реакций в виде - k*(T)^n*exp(-Ea/T), k-[см^3/с], Ea-[K]
		{
			Kf = Kchi[j][0]*pow(Tel,Kchi[j][1])*exp(-Kchi[j][2]/Tel);

			Kch[j]=Kf;
		}
		if(!strcmp(Rtype[j],"Tgas") || !strcmp(Rtype[j],"Radiat"))//запись констант реакций в виде - k*(T)^n*exp(-Ea/T), k-[см^3/с], Ea-[K]
		{
			Kf = Kchi[j][0]*pow(Tch,Kchi[j][1])*exp(-Kchi[j][2]/Tch);

			Kch[j]=Kf;
		}
		if(!strcmp(Rtype[j],"Tgas-VT"))//запись констант реакций в виде - k*(T)^n*exp(-Ea/T^0.33), k-[см^3/с], Ea-[K]
		{
			double C1;

			C1 = pow(Tch,0.33);

			Kf = Kchi[j][0]*pow(Tch,Kchi[j][1])*exp(-Kchi[j][2]/C1);

			Kch[j]=Kf;
		}
        if(!strcmp(Rtype[j],"VV"))//запись констант реакций в виде - k*(T)^n*exp(-Ea/T^0.33), k-[см^3/с], Ea-[K]
        {
            kv10 = Kchi[j][0]*pow(Tch/300.0,Kchi[j][1])*exp(-Kchi[j][2]/Tch);

            Nv = chem_VV_VT_const(Kch,N,j,"VV",kv10,Tch);

            j += Nv-1;
        }
        if(!strcmp(Rtype[j],"VV'"))//запись констант реакций в виде - k*(T)^n*exp(-Ea/T^0.33), k-[см^3/с], Ea-[K]
        {
            kv10 = Kchi[j][0]*pow(Tch/300.0,Kchi[j][1])*exp(-Kchi[j][2]/Tch);

            Nv = chem_VV_VT_const(Kch,N,j,"VV'",kv10,Tch);

            j += Nv-1;
        }
        if(!strcmp(Rtype[j],"VT"))//запись констант реакций в виде - k*(T)^n*exp(-Ea/T^0.33), k-[см^3/с], Ea-[K]
        {
            kv10 = Kchi[j][0]*pow(Tch,Kchi[j][1])*exp(-Kchi[j][2]/Tch);

            Nv = chem_VV_VT_const(Kch,N,j,"VT",kv10,Tch);

            j += Nv-1;
        }
		if(!strcmp(Rtype[j],"Rev"))//вычисление констант обратных реакций
		{
			double dH,dS,Kb,Keq,ChMm;
			double Patm = 760*p0;

			if(xRev==0)
				gas_HCpSi_calc(Tch,N);

			xRev += 1;

			dS = 0;
			dH = 0;
			for(n=0;n<N;n++)
			{
				ChMm = ChM[n][j-1]*(Mi[n]*ma);///!!!!!Mass!!!

				dH += HCpSi[0][n]*ChMm;//[эрг]*1e-10*Na;/[кДж/моль]-для проверки
				dS += HCpSi[2][n]*ChMm;//[эрг/К]*1e-7*Na;//[Дж/К*моль]-для проверки
			}

			Sum = 0;
			for(n=0;n<N;n++)
				Sum += ChM[n][j-1];

			Keq = exp(dS/kb-dH/(kb*Tch))*pow(Patm/(kb*Tch),Sum);

			Kb = Kf/Keq;

			Kch[j] = Kb;
		}

        //printf("R#%d\tKch=%.2e[cm^-3]\n",j+1,Kch[j]);

	}

	//Logging_Kch***********************************************************
	if(tic==0.0)//(dot==Ndots)
	{
		FILE *log;
		if(tic==0)
            log = fopen("Log_Kch.txt", "w");
        else
            log = fopen("Log_Kch.txt", "a+");

        fprintf(log,"(R#)\tReaction\tt=%.2e[s]\tdE,eV\n",tic);
        for(j=0; j<Nchem; j++)
            fprintf(log,"%s\t%.2e\t%.2lf\n",RName[j],Kch[j],dEch[j]/1.6e-12);
        fprintf(log,"\n\n");

		/*else
        {
            log = fopen("Log_Kch.txt", "r+");
            //rewind(log);

            char cmt[10];
            while(strcmp(cmt,";")!=0)
                fscanf(log,"%s",&cmt);
            fprintf(log,"\tt=%.2e [s] ;",dt*nt);

            cmt[0] = '\0';
            for(j=0;j<Nchem;j++)
            {
                while(strcmp(cmt,";")!=0)
                    fscanf(log,"%s",&cmt);

                fprintf(log,"\t%.2e ;",Kch[j]);
                cmt[0] = '\0';
            }
            fprintf(log,"\n\n");
        }*/

		fclose(log);
	}

}
int chem_VV_VT_const(double *Kch,int N,int j,char vin[],double kv10,double Tch)// расчёт констант скоростей колебательных реакций
{
    double dvv,Kvv,dvt,Kvt;
    double dEv;
    int Nvib = 0;
    int n,n1,n2,quant,Vmax,Vmax1,Vmax2,v1,v2;

    gas_HCpSi_calc(Tch,N);

    n1 = ChNRe[j][0]-1;
    n2 = ChNRe[j][1]-1;

    //FILE *log;
    //log = fopen("Log_VV_VT.txt", "r+");

    if(!strcmp(vin,"VV"))
    {
        n = n1;
        quant = n2-n1;

        //Gordiets et. al_"Kineticheskie proc v gasah"_1980_Nauka
        //Billing
        //Валянский_КвантЭл_1984

        dvv=-6.85/sqrt(Tch);

        Vmax = VMax[n];
        for(v1=0;v1<Vmax-quant;v1++)//(v1=0;v1<Vmax-1;v1++)
        {
            for(v2=quant;v2<Vmax;v2++)//(v2=1;v2<Vmax;v2++)
            {
                dEv = ((HCpSi[0][n+v1]-HCpSi[0][n+v1+quant])+(HCpSi[0][n+v2]-HCpSi[0][n+v2-quant]))*Mi[n]/kb;//energy_defect[K]

                Kvv = kv10*(v1+quant)*v2*exp(dvv*abs(v1-(v2-quant)))*(1.5-0.5*exp(dvv*abs(v1-(v2-quant))));//*exp(dEv*(v1-(v2-quant))/Tch);

                if(dEv<0)//need to check!!!!!!
                    Kvv *= exp(dEv/Tch);

                Kch[j+Nvib] = Kvv;

                Nvib++;

                //fprintf(log,"R#%d\t%.2lf\n",j+Nvib-1,dEv);
            }
        }

    }
    if(!strcmp(vin,"VV'"))
    {
        n2 -= 1;//
        Vmax1 = VMax[n1];
        Vmax2 = VMax[n2];

        //in analogy with VV: MUST be validated!!!
        //Gordiets et. al_"Kineticheskie proc v gasah"_1980_Nauka
        dvv=-6.85/sqrt(Tch);

        for(v1=0;v1<Vmax1-1;v1++)
        {
            for(v2=1;v2<Vmax2;v2++)
            {
                //energy_defect[K]
                dEv = (HCpSi[0][n1+v1]*Mi[n1+v1]-HCpSi[0][n1+v1+1]*Mi[n1+v1+1])/kb;
                dEv += (HCpSi[0][n2+v2]*Mi[n2+v2]-HCpSi[0][n2+v2-1]*Mi[n2+v2-1])/kb;

                Kvv = kv10*(v1+1)*v2*exp(dvv*abs(v1-(v2-1)))*(1.5-0.5*exp(dvv*abs(v1-(v2-1))));//*exp(dEv*(v1-(v2-1))/Tch);

                if(dEv<0)//need to check!!!!!!
                    Kvv = Kvv*exp(dEv/Tch);

                Kch[j+Nvib] = Kvv;
                Kch[j+Nvib+1] = Kvv*exp(-dEv/Tch);

                Nvib += 2;
            }
        }
    }
    if(!strcmp(vin,"VT"))
    {
        n1 -= 1;
        Vmax = VMax[n1];

        //Approximation!!!

        if(!strcmp(Spec[n1],"N2(0)"))
        {

            if(!strcmp(Spec[n2],"N2(0)"))
            {
                dvt = 2.87/pow(Tch,0.333333);//Billing

                for(v1=1;v1<Vmax;v1++)
                {
                    Kvt = kv10*v1*exp((v1-1)*dvt);
                    dEv = (HCpSi[0][n1+v1]-HCpSi[0][n1+v1-1])*Mi[n1]/kb;//[K]

                    Kch[j+Nvib] = Kvt;
                    Kch[j+Nvib+1] = Kvt*exp(-dEv/Tch);

                    //fprintf(log,"R#%d\t%.2lf\n",j+Nvib,dEv);
                    //fprintf(log,"R#%d\t%.2lf\n",j+Nvib+1,-dEv);

                    Nvib += 2;
                }

            }
            if(!strcmp(Spec[n2],"N"))//Capitelli,Ferreira_"Plasma_Kinetics"_2000
            {
                double Ev, beta = 0.065;//Capitelli_2000 //beta = 0.065;//Capitelli_2006
                //int dvmax = 5;

                for(v1=1;v1<Vmax;v1++)
                {
                    Ev = (HCpSi[0][n1+v1]-HCpSi[0][n1])*Mi[n1]/kb;//[K]
                    Kvt = kv10*exp(beta*Ev/Tch);

                    /*
                    //for multi-quantum transitions
                    if(v1<=dvmax)
                        Kvt = Kvt/v1;
                    else
                        Kvt = Kvt/dvmax;
                    */

                    dEv = (HCpSi[0][n1+v1]-HCpSi[0][n1+v1-1])*Mi[n1]/kb;//[K]

                    Kch[j+Nvib] = Kvt;
                    Kch[j+Nvib+1] = Kvt*exp(-dEv/Tch);

                    Nvib += 2;
                }
            }
        }
    }

    return Nvib;
}
void chem_runge_kutta4(int I,double *Ni,int Nion,int N,double *Kch,double *Sch,int Nchem,double *Wrad,double dt,double dte,double *tg_chem,double *te_chem,double tic,int dot)//расчёт по времени методом Рунге-Кутта 4 порядка
{
	int i,n,m,s,j,l;
	double al[4] = {0,dt/2,dt/2,dt};
	double Ci[N];
	double Rch[Nchem],Rch_rk[N][4],RRg;
	double M,Res,Yi,dEi;

    for(n=0;n<N;n++)
        Ci[n] = *(Ni+n*(LEN+2));

	///**********************расчёт по явной схеме(Рунге-Кутта_4)**************************
	for(i=0;i<1;i++)
	//for(i=0;i<4;i++)
	{
		if(i>0)
        {
            for(n=0;n<N;n++)
                Ci[n] = *(Ni+n*(LEN+2))+al[i]*Rch_rk[n][i-1];
        }

		///Расчет скоростей химических реакций************************************************************
		M = 0.0;
		for(n=1;n<N;n++)
			M += Ci[n];

		for(j=0;j<Nchem;j++)
		{
			Res = 1.0;
			for(l=0;l<3;l++)
			{
				if(ChNRe[j][l]==-1)
					Yi = M;
				else if(ChNRe[j][l]==0)
					Yi = 1.0;
				else
				{
					n = ChNRe[j][l]-1;
					Yi = Ci[n];
				}
				Res *= Yi;
			}

			Rch[j] = Kch[j]*Res;
		}

        ///Расчет Радиационных потерь тепла******************************************************************
        if(i==0)
        {
            *Wrad = 0.0;
            for(s=1;s<ChRad[0];s++)
            {
                m = ChRad[s];
                *Wrad += Rch[m]*dEch[m];///[erg/cm3/s]
            }
        }

		///Коэф-ты к правой части реакций(necessary processes)************************************************
		int Num,nnR,nR=0;
        int Jr,Jc,sNR;
        double SR_prod,SR_cons;
		for(n=0;n<N;n++)
		{
		    sNR = 0;
            if((i==0) && (tic/tic_ch)>3.0  && (nnR<NR))//&& (I==1)
            {
                for(nR=0;nR<NR;nR++)
                {
                    if(!strcmp(Spec[n],Spec_R[nR]))
                    {
                        Jr=0;Jc=0;SR_prod=0;SR_cons=0;
                        sNR=1;
                        break;
                    }
                }
            }

			Num = Chem[n][Nchem];
			Rch_rk[n][i] = 0.0;
			for(j=0;j<Num;j++)
			{
				m = Chem[n][j];
				RRg = ChM[n][m]*Rch[m];
				Rch_rk[n][i] += RRg;

                ///Поиск_минимального_химического_шага*******************************************************
                if((i==0) && (RRg!=0.0) && (Ci[n]>0.0))
                {
                    if(n<=Nion)
                        te_ch = fmin(fabs(Ci[n]/RRg),te_ch);///for_charged
                    else //if(n>Nion)
                        tg_ch = fmin(fabs(Ci[n]/RRg),tg_ch);///for_neutrals
                }

                ///Вычисление долей и скоростей реакций для компонент газа***************************************
                if(sNR!=0)
                {
                    if(ChM[n][m]<0)
                    {
                        RC[nR][I][Jc] = fabs(RRg);

                        Jc += 1;
                        SR_cons += fabs(RRg);
                    }
                    else
                    {
                        RP[nR][I][Jr] = RRg;

                        Jr += 1;
                        SR_prod += RRg;
                    }

                    if(j==Num-1)
                    {
                        RP[nR][I][Jr] = SR_prod;
                        RC[nR][I][Jc] = SR_cons;

                        nnR++;
                    }
                }
            }
		}
	}

	Igch += 1;
	if(Igch==LEN)
    {
        *tg_chem = tg_ch;
        tg_ch = 1.0;
        Igch = 0;
    }

    Iech += 1;
	if(Iech==LEN)
    {
        *te_chem = te_ch;
        te_ch = 1.0;
        Iech = 0;
    }

	///обновление концентраций**********************************************************
	for(n=0;n<N;n++)
	{
        //*(Sch+n*(LEN+2) = (Rch_rk[n][0]+2*Rch_rk[n][1]+2*Rch_rk[n][2]+Rch_rk[n][3])/6;
        *(Sch+n*(LEN+2)) = Rch_rk[n][0];
    }

	///Writing_Chem-Contributions******************************************************
	if(I==LEN && (tic/tic_ch)>3.0)
    {
        chem_spec_contrib(N,Nchem,1.0,tic);
        tic_ch = tic;
    }
	///******************************************************************************

}
void chem_half_implicit(double *Ni,int N0,int N,double *Kch,int Nchem,double *Wrad,double dt,double *tg_chem,double tic,int dot)//расчёт по времени неявным методом
{
	double Y0[N],Y[N];
	double Si[N],Li[N],Ki[N],Res[N],yRes;
	double Q,A,B,C,M,Y_M,Rch;

	int n,m,j,l,L,s,NJ;

	M = 0.0;
	for(n=0;n<N;n++)
    {
        Y0[n] = *(Ni+n*(LEN+2));
        M += Y0[n];
    }

	/*
	ChNRe[j][l] - //матрица реагентов
	Chem[m][j] - //матрица номеров реакций с компонентом m
	ChM[m][n] - //матрица итога реакции с компонентом m
	*/

	for(n=N0;n<N;n++)
	{
		Li[n] = 0.0;
		Ki[n] = 0.0;
		Si[n] = 0.0;
	}

	for(j=0;j<Nchem;j++)
	{
		Y_M = 1.0;
		L = 0;
		for(l=0;l<3;l++)
		{
			if(ChNRe[j][l]==-1)
				Y_M *= M;
			if(ChNRe[j][l]>0)
			{
				m = ChNRe[j][l]-1;

				for(n=N0;n<N;n++)
				{
					if(L==0)
						Res[n] = 1.0;

					if(m!=n)
						Res[n] *= Y0[m];
				}

				L += 1;
			}
		}

		for(n=N0;n<N;n++)
		{
			Res[n] *= Kch[j]*Y_M;

			NJ = ChM[n][j];

			if(NJ==-1)
				Li[n] += Res[n];
			else if(NJ==-2)
				Ki[n] += 2.0*Res[n];
			else if(NJ>0)
				Si[n] += NJ*Res[n];
		}
	}

    for(n=N0;n<N;n++)
    {
        //Q = dt*fabs(Si[n]-(Li[n]+Ki[n]*Y0[n])*Y0[n]);
		Q = Si[n]-(Li[n]+Ki[n]*Y0[n])*Y0[n];

		tg_ch = fmin(Y0[n]/Si[n],1.0/Li[n]);
		if(Ki[n]*Y0[n]!=0)
            tg_ch = fmin(tg_ch,1.0/(Ki[n]*Y0[n]));

        //if(fabs(Y0[n]/Q)>0.5*dt)
        if(tg_ch<dt || fabs(Q)*dt>=0.05*Y0[n])
        {
            ///evential
            yRes = (Y0[n]+dt*Si[n])/(1.0+(Li[n]+Ki[n]*Y0[n])*dt);
            if(yRes<0.0 || fabs(1-yRes/Y0[n])>0.05)
            {
                yRes = chem_sor_iterator(5,1.25,Y0[n],Si[n],Li[n],Ki[n],dt,1.0e-3);
                if(yRes<0.0)
                {
                    yRes = chem_newton_iterator(10,yRes,Si[n],Li[n],Ki[n],dt,1.0e-3);
                    if(yRes<0.0 && Q>0.0)
                        yRes = Y0[n]*1.01;
                    else if(yRes<0.0 && Q<0.0)
                        yRes = Y0[n]*0.99;
                }

            }

            Y[n] = yRes;
        }
        else
			Y[n] = Y0[n]+dt*Q;

        if(Y[n]<1.e-30)
			Y[n] = 0.0;

        *(Ni+n*(LEN+2)) = Y[n];
    }

    ///Поиск_минимального_химического_шага*******************************************************
    /*if((i==0) && (RRg!=0.0) && (Ci[n]>0.0))// && (n==51))
        tg_ch = fmin(fabs(Ci[n]/RRg),tg_ch);

	Igch += 1;
	if(Igch==LEN)
    {
        *tg_chem = tg_ch;
        tg_ch = 1.0;
        Igch = 0;
    }*/

    *tg_chem = 1.0e-8;

    ///Расчет Радиационных потерь тепла**********************************************************

    *Wrad = 0.0;
    for(s=1;s<ChRad[0];s++)
    {
        j = ChRad[s];

        Rch = Kch[j];
        for(l=0;l<3;l++)
		{
			if(ChNRe[j][l]==-1)
				Y_M = M;
			else if(ChNRe[j][l]==0)
				Y_M = 1.0;
			else
			{
				n = ChNRe[j][l]-1;
				Y_M = Y[n];//Y0[n];
			}
			Rch *= Y_M;
		}

        *Wrad += Rch*dEch[m];///[erg/cm3/s]
    }
}
void chem_half_modimplicit(double *Ni,int N0,int N,double *Kch,int Nchem,double *Wrad,double dt,double *tg_chem,double tic,int dot)//расчёт по времени неявным методом
{
	double Y0[N],Y[N];
	double Si[N],Li[N],Ki[N],Res[N],yRes;
	double Q,A,B,C,M,Y_M,Rch;
	int unc,nUnc[N+1];

	int i,n,m,j,l,L,s,NJ;

	M = 0.0;
	nUnc[N] = 0;
	for(n=0;n<N;n++)
    {
        Y0[n] = *(Ni+n*(LEN+2));
        M += Y0[n];
    }

    for(n=N0;n<N;n++)
    {
        nUnc[n] = n;
        nUnc[N]++;
    }
	/*
	ChNRe[j][l] - //матрица реагентов
	Chem[m][j] - //матрица номеров реакций с компонентом m
	ChM[m][n] - //матрица итога реакции с компонентом m
	*/

	s = 0;
	do
    {
        //for(n=N0;n<N;n++)
        for(i=0;i<nUnc[N];i++)
        {
            n=nUnc[i];

            Li[n] = 0.0;
            Ki[n] = 0.0;
            Si[n] = 0.0;
        }

        for(j=0;j<Nchem;j++)
        {
            Y_M = 1.0;
            L = 0;
            for(l=0;l<3;l++)
            {
                if(ChNRe[j][l]==-1)
                    Y_M *= M;
                if(ChNRe[j][l]>0)
                {
                    m = ChNRe[j][l]-1;

                    for(n=N0;n<N;n++)
                    {
                        if(L==0)
                            Res[n] = 1.0;

                        if(m!=n)
                            Res[n] *= Y0[m];
                    }

                    L += 1;
                }
            }

            //for(n=N0;n<N;n++)
            for(i=0;i<nUnc[N];i++)
            {
                n=nUnc[i];

                NJ = ChM[n][j];
                if(NJ!=0)
                {
                    Res[n] *= Kch[j]*Y_M;

                    if(NJ==-1)
                        Li[n] += Res[n];
                    else if(NJ==-2)
                        Ki[n] += 2.0*Res[n];
                    else if(NJ>0)
                        Si[n] += NJ*Res[n];
                }
            }
        }

        //for(n=N0;n<N;n++)
        for(i=0;i<nUnc[N];i++)
        {
            n=nUnc[i];

            Q = Si[n]-(Li[n]+Ki[n]*Y0[n])*Y0[n];

            tg_ch = fmin(Y0[n]/Si[n],1.0/Li[n]);
            if(Ki[n]*Y0[n]!=0)
                tg_ch = fmin(tg_ch,1.0/(Ki[n]*Y0[n]));

            if(tg_ch<dt)
            {
                if(s<5)
                    yRes = chem_sor_iterator(1,1.25,Y0[n],Si[n],Li[n],Ki[n],dt,1.0e-3);
                else
                    yRes = chem_newton_iterator(1,Y0[n],Si[n],Li[n],Ki[n],dt,1.0e-3);

                Y[n] = yRes;
            }
            else
                Y[n] = Y0[n]+dt*Q;

            if(Y[n]<1.e-30)
                Y[n] = 0.0;

            if(fabs(1-Y[n]/Y0[n])>1.0e-3)///no_converge
            {
                nUnc[unc] = n;
                unc++;
            }

             Y0[n] = Y[n];
        }
        nUnc[N] = unc;
        unc = 0;

        s++;

    }while(nUnc[N]!=0 || s<25);

    for(n=N0;n<N;n++)
        *(Ni+n*(LEN+2)) = Y[n];

    /*///Поиск_минимального_химического_шага*******************************************************
    if((i==0) && (RRg!=0.0) && (Ci[n]>0.0) && (n==51))
        tg_ch = fmin(fabs(Ci[n]/RRg),tg_ch);

	Igch += 1;
	if(Igch==LEN)
    {
        *tg_chem = tg_ch;
        tg_ch = 1.0;
        Igch = 0;
    }
    */

    *tg_chem = 5.0e-7;

    ///Расчет Радиационных потерь тепла**********************************************************

    *Wrad = 0.0;
    for(s=1;s<ChRad[0];s++)
    {
        j = ChRad[s];

        Rch = Kch[j];
        for(l=0;l<3;l++)
		{
			if(ChNRe[j][l]==-1)
				Y_M = M;
			else if(ChNRe[j][l]==0)
				Y_M = 1.0;
			else
			{
				n = ChNRe[j][l]-1;
				Y_M = Y[n];//Y0[n];
			}
			Rch *= Y_M;
		}

        *Wrad += Rch*dEch[m];///[erg/cm3/s]
    }
}
void chem_spec_contrib(int N,int Nchem,double Rth,double tic)//вывод скоростей реакций для выбранных компонент
{
	FILE *rrp;
	FILE *rcp;

	char str1[20];

	int nR,i,j,n,m,k,s,Jr,Jc,Num;
	double Rr[LEN+2],Rav;
	double Rp_av[Nchem],Rc_av[Nchem];

	for(nR=0;nR<NR;nR++)
	{
	    ///Timeline_of_Rates***********************************************************************
		sprintf(str1, "%s%s", Spec_R[nR], "_R(t)+.txt");
		rrp = fopen(str1, "a+");
		sprintf(str1, "%s%s", Spec_R[nR], "_R(t)-.txt");
		rcp = fopen(str1, "a+");

		for(n=0;n<N;n++)
		{
			if(!strcmp(Spec[n],Spec_R[nR]))
			{
				Jr = Chem[n][Nchem+1];//число реакций(приход), в которых участвует компонент
				Jc = Chem[n][Nchem+2];//число реакций(расход), в которых участвует компонент

				break;
			}
		}

        fprintf(rrp,"%.2e\t",tic);
		for(j=0;j<=Jr;j++)
			fprintf(rrp,"%.2e\t",RP[nR][1][j]);
		fprintf(rrp,"\n");

		fprintf(rcp,"%.2e\t",tic);
		for(j=0;j<=Jc;j++)
			fprintf(rcp,"%.2e\t",RC[nR][1][j]);
		fprintf(rcp,"\n");

		fclose(rrp);
		fclose(rcp);

//		double Rp_av[Jr+1],Rc_av[Jc+1];

        ///Rate_Fraction_over_length****************************************************************
        for(j=0;j<=Jr;j++)///Production
        {
            for(i=1;i<=LEN;i++)
                Rr[i] = RP[nR][i][j];
            Rav = LENaverage(Rr,Geom);

            for(i=1;i<=LEN;i++)
                Rr[i] = RP[nR][i][Jr];
            Rav *= 100/LENaverage(Rr,Geom);
            Rp_av[j] = Rav;///in_percent_%
        }

        for(j=0;j<=Jc;j++)///Consumption
        {
            for(i=1;i<=LEN;i++)
                Rr[i] = RC[nR][i][j];
            Rav = LENaverage(Rr,Geom);

            for(i=1;i<=LEN;i++)
                Rr[i] = RC[nR][i][Jc];
            Rav *= 100/LENaverage(Rr,Geom);
            Rc_av[j] = Rav;///in_percent_%
        }
        ///*****************************************************************************************

		///Spatial_Rates_Distribution***************************************************************
        sprintf(str1, "%s%s", Spec_R[nR], "_RLen+.txt");
        rrp = fopen(str1, "w");

        sprintf(str1, "%s%s", Spec_R[nR], "_RLen-.txt");
        rcp = fopen(str1, "w");

        fprintf(rrp,"r_cm\t");
        fprintf(rcp,"r_cm\t");

        Num = Chem[n][Nchem];
        k=0;
        s=0;
        for(j=0;j<Num;j++)
        {
            m = Chem[n][j];

            if(ChM[n][m]>0)
            {
                if(Rp_av[k]>=Rth)
                    fprintf(rrp,"%s_\t",RName[m]);
                k++;
            }
            else
            {
                if(Rc_av[s]>=Rth)
                    fprintf(rcp,"%s_\t",RName[m]);
                s++;
            }
        }
        fprintf(rrp,"%s_RSUM+_\n",Spec_R[nR]);
        fprintf(rcp,"%s_RSUM-_\n",Spec_R[nR]);

        fprintf(rrp,"<R,%%>\t");
        for(j=0;j<=Jr;j++)
        {
            if(Rp_av[j]>=Rth)
                fprintf(rrp,"%.2lf\t",Rp_av[j]);
        }
        fprintf(rrp,"\n");

        fprintf(rcp,"<R,%%>\t");
        for(j=0;j<=Jc;j++)
        {
            if(Rc_av[j]>=Rth)
                fprintf(rcp,"%.2lf\t",Rc_av[j]);
        }
        fprintf(rcp,"\n");

        for(i=1;i<=LEN;i++)
        {
            fprintf(rrp,"%.3lf\t",0.5*(l[i+1]+l[i]));
            for(j=0;j<=Jr;j++)
            {
                if(Rp_av[j]>=Rth)
                    fprintf(rrp,"%.2e\t",RP[nR][i][j]);
            }
            fprintf(rrp,"\n");

            fprintf(rcp,"%.3lf\t",0.5*(l[i+1]+l[i]));
            for(j=0;j<=Jc;j++)
            {
                if(Rc_av[j]>=Rth)
                    fprintf(rcp,"%.2e\t",RC[nR][i][j]);
            }
            fprintf(rcp,"\n");
        }

        fclose(rrp);
        fclose(rcp);
	}
}
double chem_sor_iterator(int I,double w,double y0,double S,double L,double K,double dt,double eps)
{
    double y,A,B,Norm;
    int i;

    i = 0;
    do
    {
        A = 1.0+dt*(L+K*y0);///impl-part
        B = y0+dt*S;///expl-part

        y = y0*(1-w) + w*B/A;

        Norm = fabs(1-y/y0);

        y0 = y;

        i++;

    }while(Norm>eps || i<I);

    if(Norm>eps)
        return -1.0;
    else
        return y;
}
double chem_newton_iterator(int I,double y0,double S,double L,double K,double dt,double eps)
{
    double y,A,B,Norm,f,F;
    int i;

    i = 0;
    do
    {
        A = 1.0+dt*(L+K*y0);///impl-part
        B = y0+dt*S;///expl-part

        f = (2.0*K*y0+L)*dt;

        F = A*y0-B;

        y = y0 - F/f;

        Norm = fabs(1-y/y0);

        y0 = y;

        i++;

    }while(Norm>eps || i<I);

    if(Norm>eps)
        return -1.0;
    else
        return y;
}
<file_sep>/Gas_calc.h
//Защита от повторного включения заголовка
#ifndef GAS_CALC_H
#define GAS_CALC_H
void gas_HCpSi_calc(double,int);
void gas_TP_calc(double *,double *,int,double *,double *,double *,double *,double *,double *,double *,double *,double *);
void gas_LenPrint(double *,double *,int,int,int,double *,double *,double *,double *,double *,double *,double *,double *,double*,double*,double*,double);
void gas_TimePrint(double *,int,double,double,double,double,double,double,double,double,double,double,double,double,double,double,double,double,double,double);
void gas_SavePrint(double *,int,double *,double *,double *,double *,double *,double);
void gas_LenAverage(int,char *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *);

#endif
<file_sep>/Chemistry_calc.h
//Защита от повторного включения заголовка
#ifndef CHEMISTRY_CALC_H
#define CHEMISTRY_CALC_H

int chem_make_react(int,char *);
void chem_read_react(int,int);
void chem_const(double *,int,int,int,double,double,double);
void chem_runge_kutta4(int,double *,int,int,double *,double *,int,double *,double,double,double *,double *,double,int);
void chem_half_implicit(double *,int,int,double *,int,double *,double,double *,double,int);
void chem_half_modimplicit(double *,int,int,double *,int,double *,double,double *,double,int);
double chem_sor_iterator(int,double,double,double,double,double,double,double);
double chem_newton_iterator(int,double,double,double,double,double,double);
void chem_spec_contrib(int,int,double,double);
int chem_VV_VT_make(int,char *,char *,char *,double *,int);
int chem_VV_VT_const(double *,int,int,char *,double,double);

#endif
<file_sep>/1DPoisson.cpp
# include "1D_MainFun.h"
# include "1DPoisson.h"
# include "1Dmesh.h"
# include "Initials.h"

extern double GF_C[LEN+2],GF_L[LEN+2],GF_R[LEN+2];
extern int vL,vR;
extern double VolL,VolR,WextL,WextR,EpsL,EpsR,DepL,DepR;

double Poisson_SORsolve(double *Fi,double *E,double Vext,double *Xi,double *Ni,double *Ngas,int Npos,int Nneg,double tic)
{
	/*
	**************************************************************

	Poisson equation solution with
	Successive over Relaxation (SOR) Method

	if w = 1.0
	SOR = Gauss-Seidel

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

    /*
	!!!!!!!!!Необходимые доработки:
        - Поток зарядов на стенку
        - улучшить алгоритм сходимости (не гонять весь диапазон по координате)_ для "частых" сеток
        - определить алгоритм выбора константы w
	*/

	/*
                left wall                                                               right wall
    Ni[i]         |                                                                       |
    Fi[i]     [0] | [1]   [2]	[3]		   		 [i-1]  [i]  [i+1]					  [I] |[I+1]
            |--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|---------->Len
    l[i]   [0]   [1]   [2]   [3]         	  [i-1]  [i]  [i+1] [i+2]              [I]  [I+1] [I+2]
    E[i]          |                                                                       |
                  |                                                                       |
    */

    int i,n;
	double Res,Ch,RHS[LEN+2];
    int cnt,Conv;//conv[LEN+2];
    double VL,VR,SigL,SigR,Norm;
    double E0[LEN+1],Fi0[LEN+2],dEf,Emax;

    ///Initial_field*********************************
    for(i=0;i<=LEN;i++)
    {
        Fi0[i] = Fi[i];
        E0[i] = E[i];
    }
    Fi0[LEN+1] = Fi[LEN+1];

    ///Volume_Charge_calculation*********************
    for(i=0;i<=LEN+1;i++)
    {
        Ch = -Ni[i] ;//electrons//Ni[i];
        for(n=1;n<=Npos;n++)//positive ions
            Ch += Ni[n*(LEN+2)+i];//Ni[n*(LEN+2)+i];

        for(n>Npos;n<=Npos+Nneg;n++)//negative ions
            Ch += -Ni[n*(LEN+2)+i] ;

        //if(fabs(Ch/Xi[i])<1.e-3)
            //Ch = 0.0;

        RHS[i] = -4*pi*e*Ch;//Ngas[i];
    }

    ///Surface_charge_density************************
    SigL = 0.0;
    SigR = 0.0;///dSig/dt = - Jrad;

    if(Vext!=0.0)///if_ExtVoltage_changes
    {
        VL = VolL;
        if(VolL!=0.0)
            VL = Vext;

        VR = VolR;
        if(VolR!=0.0)
            VR = Vext;
    }
    else///if_ExtVoltage_const
    {
        VL = VolL;
        VR = VolR;
    }
    if(WextL!=0.0)///[MHz]
        VL *= sin(2*pi*WextL*tic);
    if(WextR!=0.0)
        VR *= sin(2*pi*WextR*tic);

    ///SOR_cycle_Fi_calculation**********************
    /*cnt = 0;
    double w = 1.25;//Gauss-Seidel if w = 1.0
	do
	{
        Poisson_boundary(0,Fi,vL,VolL,EpsL,DepL,RHS[0],SigL);///Left_boundary

        Conv = 0;
        for(i=1;i<=LEN;i++)
        {
            //Res = Fi[i]*(1-w) + w*(RHS[i]-GF_L[i]*Fi[i-1]-GF_R[i]*Fi[i+1])/(-GF_C[i]);
            Res = Fi[i] + w*((RHS[i]-GF_L[i]*Fi[i-1]-GF_R[i]*Fi[i+1])/(-GF_C[i])-Fi[i]);
            ///Res = Fi[i] + w*((RHS[i]-Fi[i-1]-Fi[i+1])/(-2)-Fi[i]);///test


            //Norm = (Fi[i]-Res)*(Fi[i]-Res);
            //Norm = sqrt(Norm);
            //Norm = fabs(Res-Fi[i]);
            Norm = fabs(1-Res/Fi[i]);

            if(Norm<1.0e-3)//(fabs(Fi[i]-Res)<exact);//1.0-Res/Fi[i])<1.0e-5_i=1 &&
                Conv++;

            Fi[i] = Res;///[СГС]

            //if(Fi[i]<1.0e-5)
                //Fi[i] = 0.0;
        }

        Poisson_boundary(LEN+1,Fi,vR,VolR,EpsR,DepR,RHS[LEN+1],SigR);///Right_boundary

        cnt++;

	}while(Conv<LEN || cnt<30);//(Conv<LEN
    */

    ///SWEEP_method***************************************************************
    double  A,B,C,F,den,
            al[LEN+1],bet[LEN+1];

    //Poisson_boundary(0,Fi,vL,VL,EpsL,DepL,RHS[0],SigL);///Left_boundary
    Poisson_boundary(LEN+1,Fi,vR,VR,EpsR,DepR,RHS[LEN+1],SigR);///Right_boundary
    for(i=1;i<=LEN;i++)
    {
        ///SWEEP Coefficients:
        A = GF_L[i];///[i-1](left_cell)
        B = GF_R[i];///[i+1](right_cell)
        C = -GF_C[i];///[i](center_cell)
        F = RHS[i];///RHS-part

        if(i==1)//see SpecTransport();
        {
            //al[i-1] = 1.0;///see_SpecTransportBoundary();
            al[i-1] = 0.0;///see_SpecTransportBoundary();
            bet[i-1] = 0.0;//0.0;
        }

        den = 1.0/(A*al[i-1]+C);
        al[i] = -B*den;
        bet[i] = (F-A*bet[i-1])*den;
    }

    ///Reverse_sweep_cycle************************************************
    for(i=LEN;i>=0;i--)
        Fi[i] = al[i]*Fi[i+1]+bet[i];

    ///Potential_Smoothing???
    for(i=0;i<=LEN+1;i++)
        //Fi[i] = 0.9*Fi0[i]+0.1*Fi[i];
        Fi[i] = 0.7*Fi0[i]+0.3*Fi[i];

    ///Field_Calculation**************************************************
    for(i=0;i<=LEN;i++)
        E[i] = -2.0*(Fi[i+1]-Fi[i])/(l[i+2]-l[i]);///[СГС]

    //dEf = 0.0;
    Emax = 0.0;
    for(i=0;i<=LEN;i++)
        //dEf = fmax(fabs(E[i]-E0[i]),dEf);
        Emax = fmax(fabs(E[i]),Emax);

    return Emax;
}
void Poisson_boundary(int i,double *Fi,int v,double Vext,double Eps,double Dep,double Rhs,double Sig)
{
    int iB,iV,iV2;
    double dl,l0,l1,l2,dl2;

    if(i==0)///Left_boundary
    {
        iB = 0;///boundary
        iV = 1;///volume
        iV2 = 2;///volume2

        dl = 0.5*(l[2]-l[0]);
        l0 = 0.5*(l[1]+l[0]);
        l1 = 0.5*(l[2]+l[1]);
        l2 = 0.5*(l[3]+l[2]);
        dl2 = 0.5*(l[3]-l[1]);
    }
    else///Right_boundary
    {
        iV = LEN;///volume
        iV2 = LEN-1;///volume2
        iB = LEN+1;///volume

        dl = 0.5*(l[LEN+2]-l[LEN]);
        l0 = 0.5*(l[LEN+2]+l[LEN+1]);
        l1 = 0.5*(l[LEN+1]+l[LEN]);
        l2 = 0.5*(l[LEN]+l[LEN-1]);
        dl2 = 0.5*(l[LEN+1]-l[LEN-1]);
    }

    ///Cases:
    if(v==0)///defined
        Fi[iB] = Vext;
    if(v==1)///symmetry
        //Fi[iB] = Fi[iV];
        Fi[iB] = Fi[iV]-Rhs*l0*(l1-l0);
        //Fi[iB] = 0.99999*Fi[iV];
        //Fi[iB] = ((Fi[iV2]-Fi[iV])*l0+Fi[iV]*l2-Fi[iV2]*l1)/dl2;//linear_extrapolation

    if(v==2)///dielectric
        Fi[iB] = (Fi[iV]+Vext*Eps*dl/Dep+4.0*pi*Sig*dl)/(1.0+Eps*dl/Dep);
}

<file_sep>/1D_MainFun.h
# include <stdio.h>
# include <math.h>
# include <string.h>
# include <stdlib.h>
//# include <algorithm>

# include "Initials.h"
# include "1Dmesh.h"
# include "1DPoisson.h"
# include "1DTransport.h"
# include "Chemistry_calc.h"
# include "EEDF_calc.h"
# include "Gas_calc.h"

#ifndef MAINFUN_H
#define MAINFUN_H

# define LEN 200  //dots per axis
# define Nmax 20
# define CSmax 100
# define NRmax 200

const double
	pi = 3.141592653589,
	c = 2.997924562e+10,
	ma = 1.67e-24,//[г]
	me = 9.1e-28,//[г]
	e = 4.8e-10,//[СГС]
	eKl = 1.60217662e-19,//[Кл]
	Eabs = 300,//E[В/см]=300E[абс.СГС]
	Na = 6.022e+23,
	kb = 1.38e-16,
	eV_K = 11605,//1 эВ в кельвинах К
	p0 = 1333.22, // коэф-т перевода Торр --> СГС [эрг/см^3]
	exact = 1.0e-3,//точность счета
    cm_eV = 8065.5447;//коэф-т перевода cm-1 --> эВ


extern int NR,Ndots;//N,Nt,Nte,Nchem,Nedf,Ndots,
extern char Spec[Nmax][10],Spec_R[Nmax][10];
extern double Mi[Nmax],HCpSi[3][Nmax],CXi[Nmax][2][8];
extern double Emax,dE,dEev;
extern char Geom[10];
extern double Len,l[LEN+3],Hght;
extern double Gamma[Nmax][2],Tw;
extern double CDi[Nmax][25],CMui[Nmax];
extern char CSFile[50],ChemFile[50],EedfTblFl[50];
extern char RName[NRmax][100];

bool converge(double *,double *,double);
int sign(double);
double LENintegral(double *,char *);
double LENaverage(double *,char *);

#endif
<file_sep>/1DTransport.h
#ifndef TRANSPORT_H
#define TRANSPORT_H

void Trasport_coefs_calc(int,int,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *);
double* Transport_SWEEPsolve(double *,double *,double *,double *,double *,double *,double);
double* Transport_SWEEPsolve_mod(double *,double *,double *,double *,double *,double *,double *,double);
void SpecTransport(int,double *,double *,double *,double *,double *,double *,double *,double);
void SpecTransport_expl(double *,double *,double *,double *,double *,double *,double *,double);
void HeatTransport(double *,double *,double *,double *,double *,double *,double *,int,double *,double *,double *,double);
void TeTransport(double *,double *,double *,double *,int,double);
void TeTransport_expl(double *,double *,double *,double *,double *,double *,double *,double *,double);
void TransportBoundary(int,double *,double *,double *,double *,double *,double *,double *,double);
void TransportBoundary_mod(int,int,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double);

#endif // TRANSPORT_H
<file_sep>/1Dmesh.h
#ifndef MESH_H
#define MESH_H

void Mesh_gen(double);
void Mesh_GFcalc(char *);

#endif // MESH_H
<file_sep>/EEDF_calc.h
//Защита от повторного включения заголовка
#ifndef EEDF_CALC_H
#define EEDF_CALC_H
int EEDF_read_CS(int,char *);
void EEDF_calc(int,double *,double *,double *,double,double,double,double *,double *,double *,double *,double,double,int,double);//решение уравнения Больцмана
void EEDF_const_calc(char *,int,double *,double,double,double,double *,double *,double *,double *,double);
void EEDF_print(double *,double,double,double,double,char *);
void EEDF_table_calc(int,int,int,double *,double,double,double,char *);
void EEDF_table_print(int,int,char *);
int EEDF_table_read(int,int,char *);
#endif
<file_sep>/1DTransport.cpp
# include "1D_MainFun.h"
# include "1DTransport.h"
# include "Gas_calc.h"
# include "1Dmesh.h"

extern double GF_C[LEN+2],GF_L[LEN+2],GF_R[LEN+2];
double al_bound[Nmax+1][2],bet_bound[Nmax+1][2];

void Trasport_coefs_calc(int N,int Nion,double *Ni,double *E,double *Del,double *Muel,double *Di,double *Mui,double *Lam,double *Pgas,double *Tgas,double *Ngas,double *Te)
{
    double Damb,SumNi,DPT,MuPTi,EN;
    double A,B,xDiff,Dix,DN;
    int i,n,ni,k,ki,s,S;
    int Z[Nion+1];

    //Tgas[0] = Tgas[1];
    Te[0] = Te[1];
    Te[LEN+1] = Te[LEN];

    Z[0] = -1;
    for(n=1;n<=Nion;n++)
        Z[n] = sign(CMui[n]);

    for(i=0;i<=LEN+1;i++)
    {
        //Tgas[i] = 415.9;
        //Ngas[i] = 3.22e16*300/Tgas[i];
        EN = fabs(E[i])*1.0e17*Eabs/Ngas[i];///[Td]//5;//91.52;//

        for(n=0;n<N;n++)
        {
            ni = n*(LEN+2)+i;

            if(n==0)
            {
                Di[ni] = Del[i];
                Mui[ni] = Z[n]*Muel[i];//0.0;//

                *(Muel+ni+(LEN+2)) *= Z[n];///MuEel
                if(i==LEN+1)
                {
                    Di[0] = Di[1];
                    *(Del+(LEN+2)) = *(Del+(LEN+2)+1);
                    Di[LEN+1] = Di[LEN];
                    *(Del+2*(LEN+1)+1) = *(Del+2*(LEN+1));//*(Del+(LEN+2)+(LEN+1)) = *(Del+(LEN+2)+LEN);

                    Mui[0] = Mui[1];
                    *(Muel+(LEN+2)) = *(Muel+(LEN+2)+1);
                    Mui[LEN+1] = Mui[LEN];
                    *(Muel+2*(LEN+1)+1) = *(Muel+2*(LEN+1));
                }
            }
            else if(n>0 && n<=Nion)
            {
                DPT = pow(Tgas[i],CDi[n][1])/Pgas[i];
                MuPTi = 1.0;///

                ///Di(E/N)_calculation**********************************************
                Dix = CDi[n][0]*DPT;
                //DN = Dix*Ngas[i];
                if(CDi[n][3]!=0.0 && EN>=0.0)///Di(Ez/N)_dependence
                {
                    S = CDi[n][2]-1;

                    if(EN>CDi[n][S*2+3])
                        Dix  *= CDi[n][S*2+4];
                    else
                    {
                        s=0;
                        while(EN>CDi[n][s*2+3] && s<=S)
                            s++;

                        if(s==0)
                        {
                            A = (CDi[n][s*2+4]-1.0)/CDi[n][s*2+3];
                            B = 1.0;
                        }
                        else
                        {
                            A = (CDi[n][s*2+4]-CDi[n][(s-1)*2+4])/(CDi[n][s*2+3]-CDi[n][(s-1)*2+3]);
                            B = CDi[n][(s-1)*2+4]-A*CDi[n][(s-1)*2+3];
                        }
                        xDiff = A*EN+B;
                        Dix  *= xDiff;
                    }
                }
                Di[ni] = Dix;

                ///Mui(E/N)_calculation****************************************
                Mui[ni] = CMui[n]*MuPTi;
            }
            else
            {
                DPT = pow(Tgas[i],CDi[n][1])/Pgas[i];
                Di[ni] = CDi[n][0]*DPT;
                Mui[ni] = 0.0;//1.0;//e/me/Vm_av;
            }

            ///Lam = 1.0e+7*(2.714+5.897e-3*Tin)*1.0e-4; //[Эрг/(с*см*K)]Ref. - Nuttall_JRNBS_1957 (N2)
            if(n==N-1)
                //Lam[i] = 34*pow(Tgas[i],0.76);///[erg/cm*s*K]//[N2-Mankelevich_form_consists_with_Справочник-ФизВеличин]
                Lam[i] = 33*pow(Tgas[i],0.78);///[erg/cm*s*K]//[O2-My_approx_from_data_Справочник-ФизВеличин]
                //Lam[i] = 2.4/100*(6.025e-5*Tgas[i]+8.55e-3)*1e7;///[erg/cm*s*K]//[O2-Proshina_data]_=My*2.3
                //Lam[i] = 32*pow(Tgas[i],0.71);///[erg/cm*s*K]//[Ar-My_approx_from_data_Справочник-ФизВеличин]
        }
    }
}
void SpecTransport(int n,double *Xi,double *Ni,double *Di,double *Mui,double *Ngas,double *E,double *Si,double dt)
{
	/*
                left wall                                                               right wall
    Ni[i]         |                                                                       |
    Fi[i]     [0] | [1]   [2]	[3]		   		 [i-1]  [i]  [i+1]					  [I] |[I+1]
            |--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|---------->Len
    l[i]   [0]   [1]   [2]   [3]         	  [i-1]  [i]  [i+1] [i+2]              [I]  [I+1] [I+2]
    E[i]          |                                                                       |
                  |                                                                       |
    */

    double Vdr[LEN+1];
    double *res;

    ///**************Defining_Drift-coefficient and Source_term******************
    int i;
    for(i=1;i<=LEN;i++)
    {
        if(i==1)
        {
            Vdr[i-1] = 0.5*(Mui[i]+Mui[i-1])*E[i-1];

            Si[i-1] = 0.0;
            Si[LEN+1] = 0.0;
        }

        Vdr[i] = 0.5*(Mui[i+1]+Mui[i])*E[i];

        //Si[i] = 0.0;
    }

    ///New_Jd=-D*N*dXi/dx_Vdr=N*Xi*********************************
    res = Transport_SWEEPsolve_mod(Ngas,Xi,Di,Vdr,Si,&al_bound[n][0],&bet_bound[n][0],dt);

    for(i=1;i<=LEN;i++)
    {
        Xi[i] = *(res+i);
        Ni[i] = *(res+i)*Ngas[i];///reNEW

        /*if(Ni[i]<1.0e-30)
            Ni[i] = 1.e-30;*/
    }
}
double* Transport_SWEEPsolve_mod(double *Ng,double *Xi,double *Di,double *Vi,double *Si,double *al_b,double *bet_b,double dt)
{
    /*
    //-----------------------------------------------------------------------------------------------------------------------------

	Transport equation (in Drift-Diffusion approximation ) solution with implicit SWEEP/Shuttle Method.

	dNi/dt = d(D*(dNi/dx))/dx - d(Vi*Ni)dx + S , S - source term.

	Drift-Diffusion Term (DDT) is implemented using exact solution for steady 1D DDT problem:

	d(RoV*F)/dx = dd(D*F)/dx*dx         (*),

	(F(x)-F(0))/(F(L)-F(0))=(exp(Pe*x/L)-1)/(exp(Pe)-1), Pe=RoV*L/D - Peclet Number

    In our Case:

    1)Species:
    J[i+1/2] = Vdr[i+0.5]*(N[i]-(N[i+1]-N[i])/(exp(Pe[i+0.5])-1));

    "Power Law Scheme" (Pantakar_1980_page-95) used for exp() approximation

    Linearization of Equation (*):

    J[i+1/2] = Vdr[i+0.5]*(N[i]-(1-0.1*Pe)^5(N[i+1]-N[i])/Pe), for 0<Pe<=5;

       [i-1]  [i]  [i+1]
    --|--x--|--x--|--x--|--

    A[i]*N[i] = A[i+1]*N[i+1] + A[i-1]*N[i-1]

    Coefficients (linear geometry):
    A[i+1] = D[i+0.5]*F(|Pe[i+0.5]|)+Max(-Vdr[i+0.5];0)
    A[i-1] = D[i-0.5]*F(|Pe[i-0.5]|)+Max(Vdr[i-0.5];0)
    A[i] = A[i+1]+A[i-1]+(Vdr[i+0.5]-Vdr[i-0.5])

    where:
    F(|Pe|) = Max(0;(1-0.1*|Pe|)^5)

    In "RADIAL" case:

    AA[i+1] = A[i+1]*r[i+0.5]/(r[i]*(r[i+0.5]-r[i-0.5]))
    AA[i-1] = A[i-1]*r[i-0.5]/(r[i]*(r[i+0.5]-r[i-0.5]))
    AA[i] = AA[i+1]+AA[i-1]+(r[i+0.5]*Vdr[i+0.5]-r[i-0.5]*Vdr[i-0.5])/(r[i]*(r[i+0.5]-r[i-0.5]))

    Chemistry part considered as explicit (source term)

	//-----------------------------------------------------------------------------------------------------------------------------
	*/

	double  D_L,D_R,V_L,V_R,F_L,F_R,dL,
            Pe,fPe,
            A_L,A_R,A_C;

    double  A,B,C,F,den,
            al[LEN+1],bet[LEN+1];

    ///SWEEP-SHUTTLE_CICLE*****************************************

    ///Defining_sweep_coefficients
    int i;
    for(i=1;i<=LEN;i++)
    {
        ///Defining_sweep_coefficients
        ///Diffusion
        D_L = D_R;
        if(i==1)
            D_L = 0.5*(Ng[i]+Ng[i-1])*0.5*(Di[i]+Di[i-1]);//0.5*(Di[i]*Ng[i]+Di[i-1]*Ng[i-1]);//
        D_R = 0.5*(Ng[i]+Ng[i+1])*0.5*(Di[i]+Di[i+1]);//0.5*(Di[i+1]*Ng[i+1]+Di[i]*Ng[i]);//

        ///Drift
        V_L = V_R;
        if(i==1)
            V_L = Vi[i-1]*0.5*(Ng[i]+Ng[i-1]);//0.5*(Vi[i]+Vi[i-1]);
        V_R = Vi[i]*0.5*(Ng[i+1]+Ng[i]);//0.5*(Vi[i+1]+Vi[i]);

        /*//Drift_NEW????
        V_L = V_R;
        if(i==1)
        {
            if(Vi[i-1]>0.0)
                V_L = Vi[i-1]*Ng[i-1];
            else
                V_L = Vi[i-1]*Ng[i];//0.5*(Vi[i]+Vi[i-1]);
        }
        if(Vi[i]>0.0)
            V_R = Vi[i]*Ng[i];
        else
            V_R = Vi[i]*Ng[i+1];
        */

        ///DDT_Coefficients:
        ///right-edge
        dL = 0.5*(l[i+2]-l[i]);
        F_R = V_R*dL;
        Pe = fabs(F_R/D_R);
        //Power-Law_scheme:
        fPe = pow((1.0-0.1*Pe),5.0);
        A_R = D_R*fmax(0.0,fPe)+fmax(-F_R,0.0);
        //Exp-scheme:
        //A_R = V_R/(exp(Pe)-1.0);

        ///left-edge
        dL = 0.5*(l[i+1]-l[i-1]);
        F_L = V_L*dL;
        Pe = fabs(F_L/D_L);
        //Power-Law_scheme:
        fPe = pow((1.0-0.1*Pe),5.0);
        A_L = D_L*fmax(0.0,fPe)+fmax(F_L,0.0);
        //Exp-scheme:
        //A_L = V_L/(exp(Pe)-1.0);

        ///Accounting for Geometry Factors_POWER-scheme:
        A_R = A_R*GF_R[i];
        A_L = A_L*GF_L[i];
        A_C = A_R+A_L+GF_R[i]*F_R-GF_L[i]*F_L;

        /*
        ///Accounting for Geometry Factors_EXP-scheme:
        A_R = A_R*GF_R[i];
        A_L = (A_L+V_L)*GF_L[i];
        A_C = (A_R+V_R)*GF_R[i]+A_L*GF_L[i];
        */

        ///SWEEP Coefficients_NEW-dt:
        A = -A_L*dt;///[i-1](left_cell)
        B = -A_R*dt;///[i+1](right_cell)
        C = A_C*dt+Ng[i];///[i](center_cell)
        F = Xi[i]*Ng[i]+Si[i]*dt;///RHS-part

        /*///SWEEP Coefficients_NEW-Stationary_do_not_work??:
        A = -A_L;///[i-1](left_cell)
        B = -A_R;///[i+1](right_cell)
        C = A_C;///[i](center_cell)
        F = Si[i];///RHS-part*/

        if(i==1)//see SpecTransport();
        {
            al[i-1] = al_b[0];///see_SpecTransportBoundary();
            bet[i-1] = bet_b[0];//0.0;
        }

        den = 1.0/(A*al[i-1]+C);
        al[i] = -B*den;
        bet[i] = (F-A*bet[i-1])*den;

        /*if(i==LEN)//see SpecTransport();
        {
            al[i] = al_b[1];///see_SpecTransportBoundary();
            bet[i] = bet_b[1];
        }*/

    }

    //Reverse_sweep_cycle************************************************
    for(i=LEN;i>=0;i--)
    {
        Xi[i] = al[i]*Xi[i+1]+bet[i];
        if(Xi[i]*Ng[i]<1.e-30)
            Xi[i] = 0.0;
    }

    return &Xi[0];
}
void TransportBoundary_mod(int N,int Nion,double *Ni,double *Xi,double *Di,double *Mui,double *E,double *Ngas,double *Nel,double *Tgas,double *Pgas,double *Tv,double *Te,double Twall)
{
    int n,nX=N-1,nR,nL;//N-1
    double Pe,Ped,D,Vt,Vd,JL,JR;
    double al,alXL,alXR,dlL,dlR,DXL,DXR,dNL,dNR;
    double Nnew,sign;

    dlL = 0.5*(l[2]-l[0]);
    dlR = 0.5*(l[LEN+2]-l[LEN]);

    ///Temperature_Boundary_conditions*****************************************
    ///Left_boundary**********************
    if(!strcmp(Geom,"axial"))
        Tgas[0] = Tgas[1];
    else if(!strcmp(Geom,"cartesian"))
        Tgas[0] = Twall;
    Ngas[0] = Pgas[0]/(kb*Tgas[0]);

    ///Right_boundary*********************
    Tgas[LEN+1] = Twall;
    Ngas[LEN+1] = Pgas[LEN+1]/(kb*Tgas[LEN+1]);

    ///Species_Boundary_conditions*********************************************
    dNL = 0.0;
    dNR = 0.0;
    for(n=0;n<N;n++)
    {
        nL = n*(LEN+2);
        nR = nL+LEN+1;

        ///Left_boundary***********************************************
        if(n<=Nion)///Charged_particles
        {
            if(Gamma[n][0] == 0.0)
            {
                Xi[nL] = Xi[nL+1];

                al_bound[n][0] = 1.0;
                bet_bound[n][0] = 0.0;
            }
            else if(Gamma[n][0] == 1.0)
            {
                Xi[nL] = 1.e-50;

                al_bound[n][0] = 0.0;
                bet_bound[n][0] = 0.0;
            }
            else
            {
                D = 0.5*(Di[nL+1]*Ngas[1]+Di[nL]*Ngas[0]);
                Vd = 0.5*(Mui[nL+1]*Ngas[1]+Mui[nL]*Ngas[0])*E[0];
                JL = -D*(Xi[nL+1]-Xi[nL])+Vd*Xi[nL+1];

                if(Xi[nL]-Xi[nL+1]>0.0 || JL>0.0)
                {
                    //Xi[nL] = Xi[nL+1];

                    //al_bound[n][0] = 1.0;
                }
                else
                {
                    Ped = dlR*fabs(Vd)/D;

                    if(n==0)
                        Vt = sqrt(8*kb*Te[0]*eV_K/(pi*Mi[n]));
                    else
                        Vt = sqrt(8*kb*Tgas[0]/(pi*Mi[n]));
                    Vt  *= Ngas[0];
                    Pe = dlL*Vt/D;

                    al = (1.0+sign*0.25*Gamma[n][0]*Pe)/(1.0+Ped);

                    if(n==nX)
                        alXL = al;

                    Nnew = Xi[nL+1]/al;

                    if(n!=0)
                        dNL += 0.25*Gamma[n][0]*Vt*Xi[nL];//fabs(Nnew-Nold);//

                    Xi[nL] = Nnew;

                    al_bound[n][0] = 1.0/al;
                }

                bet_bound[n][0] = 0.0;
            }
        }
        else///Neutral_particles
        {
            if(Gamma[n][0] == 0.0)//Ni[n][0] = Ni[n][1];
            {
                Xi[nL] = Xi[nL+1];

                al_bound[n][0] = 1.0;
                bet_bound[n][0] = 0.0;

                if(n==nX)
                {
                    DXL = 0.5*(Di[nL+1]*Ngas[1]+Di[nL]*Ngas[0]);
                    alXL = 1.0;
                }
            }
            else///equation for DDT with kinetic wall flux
            {
                if(Xi[nL+1]-Xi[nL]>0.0)
                    sign = 1.0;
                else
                    sign = 0.0;

                D = 0.5*(Di[nL+1]*Ngas[1]+Di[nL]*Ngas[0]);
                if(n==nX)
                    DXL = D;

                Vt = sqrt(8*kb*Tgas[0]/(pi*Mi[n]))*Ngas[0];
                Pe = dlL*Vt/D;

                al = (1.0+sign*0.25*Gamma[n][0]*Pe);

                if(n==nX)
                    alXL = al;

                Nnew = Xi[nL+1]/al;

                if(n!=nX)
                    dNL += 0.25*Gamma[n][0]*Vt*Xi[nL];//fabs(Nnew-Nold);//

                Xi[nL] = Nnew;

                al_bound[n][0] = 1.0/al;
                bet_bound[n][0] = 0.0;
            }
        }

        ///Right_boundary************************************************
        if(n<=Nion)///Charged_particles
        {
            if(Gamma[n][1] == 0.0)
            {
                Xi[nR] = Xi[nR-1];

                al_bound[n][1] = 1.0;
                bet_bound[n][1] = 0.0;
            }
            else if(Gamma[n][1] == 1.0)
            {
                Xi[nR] = 1.e-30;

                al_bound[n][1] = 0.0;
                bet_bound[n][1] = 0.0;
            }
            else
            {
                D = 0.5*(Di[nR]*Ngas[LEN+1]+Di[nR-1]*Ngas[LEN]);
                Vd = 0.5*(Mui[nR]*Ngas[LEN+1]+Mui[nR-1]*Ngas[LEN])*E[LEN];
                JR = -D*(Xi[nR]-Xi[nR-1])+Vd*Xi[nR-1];

                if(Xi[nR]-Xi[nR-1]>0.0 || JR<0.0)
                {
                    Xi[nR] = Xi[nR-1];

                    al_bound[n][1] = 1.0;
                }
                else
                {
                    Ped = dlR*Vd/D;

                    if(n==0)
                        Vt = sqrt(8*kb*Te[LEN+1]*eV_K/(pi*Mi[n]));
                    else
                        Vt = sqrt(8*kb*Tgas[LEN+1]/(pi*Mi[n]));
                    Vt  *=  Ngas[LEN+1];
                    Pe = dlR*Vt/D;

                    al = (1.0+0.25*Gamma[n][1]*Pe)/(1.0+Ped);

                    Nnew = Xi[nR-1]/al;

                    if(n!=0)
                        dNR += 0.25*Gamma[n][1]*Vt*Xi[nR];

                    Xi[nR] = Nnew;

                    al_bound[n][1] = al;
                }

                bet_bound[n][1] = 0.0;
            }
        }
        else///Neutral_particles
        {
            if(Gamma[n][1] == 0.0)//Ni[n][LEN+1] = Ni[n][LEN];
            {
                Xi[nR] = Xi[nR-1];

                al_bound[n][1] = 1.0;
                bet_bound[n][1] = 0.0;

                if(n==nX)
                {
                    DXR = 0.5*(Di[nR]*Ngas[LEN+1]+Di[nR-1]*Ngas[LEN]);

                    alXR = 1.0;
                }
            }
            else///equation for DDT with kinetic wall flux
            {
                if(Xi[nR]-Xi[nR-1]>0.0)
                    sign = 0.0;//-1.0;
                else
                    sign = 1.0;

                D = 0.5*(Di[nR]*Ngas[LEN+1]+Di[nR-1]*Ngas[LEN]);
                if(n==nX)
                    DXR = D;

                Vt = sqrt(8*kb*Tgas[LEN+1]/(pi*Mi[n]))*Ngas[LEN+1];
                Pe = dlR*Vt/D;

                al = (1.0+sign*0.25*Gamma[n][1]*Pe);

                if(n==nX)
                    alXR = al;

                Nnew = Xi[nR-1]/al;

                if(n!=nX)
                    dNR += 0.25*Gamma[n][1]*Vt*Xi[nR];

                Xi[nR] = Nnew;

                al_bound[n][1] = al;
                bet_bound[n][1] = 0.0;
            }
        }

        Ni[nL] = Xi[nL]*Ngas[0];
        Ni[nR] = Xi[nR]*Ngas[LEN+1];

        if(n==0)
        {
            Nel[nL] = Ni[nL];
            Nel[nR] = Ni[nR];
        }
    }

    ///for_Ground_state_n=nX:
    nL = nX*(LEN+2);
    nR = nL+LEN+1;
    if(dNL!=0.0)
    {
        ///OLD_version
        //bet_bound[n][0] = dNL*dlL/(Ngas[0]*Mi[n]*DXL*alXL);//

        ///NEW_version
        bet_bound[nX][0] = 0.0;//dNL*dlL/(DXL*alXL);//
        Xi[nL] = Xi[nL+1];//alXL+bet_bound[n][0];??????
        Ni[nL] = Xi[nL]*Ngas[0];
    }
    if(dNR!=0.0)
    {
        ///OLD_version
        //bet_bound[n][1] = -1.0*dNR*dlR/(Ngas[LEN+1]*Mi[n]*DXR);
        //Xi[n*(LEN+2)+LEN+1] += -1.0*bet_bound[n][1]/alXR;

        ///NEW_version
        bet_bound[nX][1] = -dNR*dlR/DXR;//
        Xi[nR] = (Xi[nR-1]-bet_bound[nX][1])/alXR;
        Ni[nR] = Xi[nR]*Ngas[LEN+1];
    }

    ///Ne-Te_Boundary_conditions***********************************************
    ///left_boundary**********************
    Te[0] = Twall/eV_K;//Te[1];
    Tv[0] = Tv[1];

    al_bound[N+1][0] = 1.0;
    bet_bound[N+1][0] = 0.0;

    ///Right_boundary**********************
    Te[LEN+1] = Twall/eV_K;//Te[LEN];
    Tv[LEN+1] = Tv[LEN];

    al_bound[N+1][1] = 1.0;
    bet_bound[N+1][1] = 0.0;

}
void HeatTransport(double *Hgas,double *Ngas,double *Tgas,double *Lam,double *Ni,double *Xi,double *Di,int N,double *E,double *J,double *Wrad,double dt)
{
    //Сетка по длине:
	/*
                left wall                                                               right wall
    Ni[i]         |                                                                       |
    Fi[i]     [0] | [1]   [2]	[3]		   		 [i-1]  [i]  [i+1]					  [I] |[I+1]
            |--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|---------->Len
    l[i]   [0]   [1]   [2]   [3]         	  [i-1]  [i]  [i+1] [i+2]              [I]  [I+1] [I+2]
    E[i]          |                                                                       |
                  |                                                                       |
    */

    //Учет нагрева электронами в упругих соударениях + нагрев газа полем + теплопроводность на стенку трубки
    //Lam = 1.0e+7*(2.714+5.897e-3*Tin)*1.0e-4; //[Эрг/(с*см^2*K)]Ref. - Nuttall_JRNBS_1957 (N2)
    //Hin += (Qel)*dt - Lam*(Tin-Tw)/pow(Rad/2.4,2.0)*dt;//w/o QE

    double  Lm,QL,QR,JHL,JHR,T0,Dn;
    int i,n,ni;

    //Calculation_of_enthalpy*************************************

    for(i=1;i<=LEN;i++)
    {
        //Учет теплопроводности******************

        QL = QR;
        if(i==1)
        {
            Lm = 0.5*(Lam[i]+Lam[i-1]);
            QL = - Lm*(Tgas[i]-Tgas[i-1]);
        }
        Lm = 0.5*(Lam[i+1]+Lam[i]);
        QR = - Lm*(Tgas[i+1]-Tgas[i]);

        //Учет диффузионного переноса************

        JHL = JHR;
        if(i==1)
        {
            T0 = (Tgas[i]+Tgas[i-1])*0.5;
            gas_HCpSi_calc(T0,N);

            JHL = 0.0;
            for(n=1;n<N;n++)//1
            {
                ni = n*(LEN+2)+i;

                Dn = 0.5*(Di[ni]*Ngas[i]+Di[ni-1]*Ngas[i-1]);
                JHL += - HCpSi[0][n]*Mi[n]*Dn*(Xi[ni]-Xi[ni-1]);
            }
        }
        JHR = 0.0;
        T0 = (Tgas[i+1]+Tgas[i])*0.5;
        gas_HCpSi_calc(T0,N);
        for(n=1;n<N;n++)//1
        {
            ni = n*(LEN+2)+i;

            Dn = 0.5*(Di[ni+1]*Ngas[i+1]+Di[ni]*Ngas[i]);
            JHR += - HCpSi[0][n]*Mi[n]*Dn*(Xi[ni+1]-Xi[ni]);
        }

        Hgas[i] += -dt*(GF_R[i]*(QR+JHR)-GF_L[i]*(QL+JHL))+dt*(fabs(J[i]*E[i])-Wrad[i]);///{J*E}=[erg/cm3/c]
    }
}
void SpecTransport_expl(double *Xi,double *Ni,double *Di,double *Mui,double *Sch,double *Ngas,double *E,double dt)
{
    int i;
    double Flux[LEN+1];
    double  D,Vdr,A,F,dL,Pe,fPe;

    ///Flux_calculations_Jd=-D*N*dXi/dx_Vdr=N*Xi*********************************
    for(i=1;i<=LEN;i++)
    {
        if(i==1)
        {
            ///DDT_Coefficients:
            //D = 0.5*(Di[1]*Ngas[1]+Di[0]*Ngas[0]);//0.5*(Ngas[1]+Ngas[0])*0.5*(Di[1]+Di[0]);//0.5*(Di[i]*Ng[i]+Di[i-1]*Ng[i-1]);//
            D = 0.5*(Ngas[1]+Ngas[0])*0.5*(Di[1]+Di[0]);

            //Vdr = 0.5*(Mui[1]*Ngas[1]+Mui[0]*Ngas[0])*E[0];
            Vdr = 0.5*(Mui[1]+Mui[0])*0.5*(Ngas[1]+Ngas[0])*E[0];

            /*dL = 0.5*(l[2]-l[0]);
            if(Vdr>=0.0)
                Flux[0] = -D*(Xi[1]-Xi[0])+Vdr*dL*Ngas[0]*Xi[0];
            else
                Flux[0] = -D*(Xi[1]-Xi[0])+Vdr*dL*Ngas[1]*Xi[1];
            */

            ///DDT_Flux:
            if(Vdr!=0.0)///DDT_Flux
            {
                dL = 0.5*(l[2]-l[0]);
                F = Vdr*dL;
                Pe = fabs(F/D);
                //Power-Law_scheme:
                fPe = pow((1.0-0.1*Pe),5.0);
                A = D*fmax(0.0,fPe)+fmax(F,0.0);

                Flux[0] = -A*(Xi[1]-Xi[0])+F*Xi[0];
            }
            else///Only_Diffusion_Flux
                Flux[0] = -D*(Xi[1]-Xi[0]);
        }

        ///DDT_Coefficients:
        //D = 0.5*(Di[i+1]*Ngas[i+1]+Di[i]*Ngas[i]);//0.5*(Ngas[i]+Ngas[i+1])*0.5*(Di[i]+Di[i+1]);
        D = 0.5*(Ngas[i]+Ngas[i+1])*0.5*(Di[i]+Di[i+1]);

        //Vdr = 0.5*(Mui[i+1]*Ngas[i+1]+Mui[i]*Ngas[i])*E[i];//0.5*(Mui[i+1]+Mui[i])*E[i];
        Vdr = 0.5*(Mui[i+1]+Mui[i])*0.5*(Ngas[i]+Ngas[i+1])*E[i];

        ///DDT_Flux:
        if(Vdr!=0.0)///DDT_Flux
        {
            dL = 0.5*(l[i+2]-l[i]);
            F = Vdr*dL;
            Pe = fabs(F/D);
            //Power-Law_scheme:
            fPe = pow((1.0-0.1*Pe),5.0);
            A = D*fmax(0.0,fPe)+fmax(-F,0.0);

            Flux[i] = -A*(Xi[i+1]-Xi[i])+F*Xi[i];
        }
        else///Only_Diffusion_Flux
            Flux[i] = -D*(Xi[i+1]-Xi[i]);
    }

    ///DDT_Equation_Solving****************************************************
    for(i=1;i<=LEN;i++)
    {
        Ni[i] += -dt*(GF_R[i]*Flux[i]-GF_L[i]*Flux[i-1])+dt*Sch[i];

        /*if(Ni[i]<1.0e-30)
            Ni[i] = 1.e-30;*/
    }
}
void TeTransport_expl(double *Te,double *Nel,double *DE,double *MuE,double *Je,double *E,double *SE,double *dTe,double dt)
{
    int i;
    double Flux[LEN+1];
    double Te0,D,Vdr,A,F,dL,Pe,fPe;
    double NE[LEN+2];

    ///NeTe_Recalculation********************************************************
    for(i=0;i<=LEN+1;i++)
        NE[i] = 3.0/2.0*Nel[i]*Te[i];///

    ///Flux_calculations_Jd=-D*N*dXi/dx_Vdr=N*Xi*********************************
    for(i=1;i<=LEN;i++)
    {
        if(i==1)
        {
            ///DDT_Coefficients:
            D = 0.5*(DE[1]+DE[0]);
            Vdr = 0.5*(MuE[1]+MuE[0])*E[0];///MuE<-->sign(-1)

            /*dL = 0.5*(l[2]-l[0]);
            if(Vdr>=0.0)
                Flux[0] = -D*(Xi[1]-Xi[0])+Vdr*dL*Ngas[0]*Xi[0];
            else
                Flux[0] = -D*(Xi[1]-Xi[0])+Vdr*dL*Ngas[1]*Xi[1];
            */

            ///DDT_Flux:
            if(Vdr!=0.0)///DDT_Flux
            {
                dL = 0.5*(l[2]-l[0]);
                F = Vdr*dL;
                Pe = fabs(F/D);
                //Power-Law_scheme:
                fPe = pow((1.0-0.1*Pe),5.0);
                A = D*fmax(0.0,fPe)+fmax(F,0.0);

                Flux[0] = -A*(NE[1]-NE[0])+F*NE[0];
            }
            else///Only_Diffusion_Flux
                Flux[0] = -D*(NE[1]-NE[0]);
        }

        ///DDT_Coefficients:
        D = 0.5*(DE[i]+DE[i+1]);
        Vdr = 0.5*(MuE[i+1]+MuE[i])*E[i];

        ///DDT_Flux:
        if(Vdr!=0.0)///DDT_Flux
        {
            dL = 0.5*(l[i+2]-l[i]);
            F = Vdr*dL;
            Pe = fabs(F/D);
            //Power-Law_scheme:
            fPe = pow((1.0-0.1*Pe),5.0);
            A = D*fmax(0.0,fPe)+fmax(-F,0.0);

            Flux[i] = -A*(NE[i+1]-NE[i])+F*NE[i];
        }
        else///Only_Diffusion_Flux
            Flux[i] = -D*(NE[i+1]-NE[i]);
    }

    ///DDT_Equation_Solving****************************************************
    for(i=1;i<=LEN;i++)
    {
        Te0 = Te[i];

        NE[i] += -dt*(GF_R[i]*Flux[i]-GF_L[i]*Flux[i-1]);//+dt*(fabs(Je[i]*E[i])/1.602e-12);//-SE[i]);///{J*E}=[erg/cm3/c]_1.602e-12;//[eV]=1.602e-12[erg]

        Te[i] = 2.0/3.0*NE[i]/Nel[i];///Te_Recalculation***********************

        dTe[i] = fabs(Te[i]-Te0);
    }
}
///******************************************NOT USED********************************************

double* Transport_SWEEPsolve(double *NNi,double *Di,double *Vi,double *Si,double *al_b,double *bet_b,double dt)
{
    /*
    //-----------------------------------------------------------------------------------------------------------------------------

	Transport equation (in Drift-Diffusion approximation ) solution with implicit SWEEP/Shuttle Method.

	dNi/dt = d(D*(dNi/dx))/dx - d(Vi*Ni)dx + S , S - source term.

	Drift-Diffusion Term (DDT) is implemented using exact solution for steady 1D DDT problem:

	d(RoV*F)/dx = dd(D*F)/dx*dx         (*),

	(F(x)-F(0))/(F(L)-F(0))=(exp(Pe*x/L)-1)/(exp(Pe)-1), Pe=RoV*L/D - Peclet Number

    In our Case:

    1)Species:
    J[i+1/2] = Vdr[i+0.5]*(N[i]-(N[i+1]-N[i])/(exp(Pe[i+0.5])-1));

    "Power Law Scheme" (Pantakar_1980) used for exp() approximation

    Linearization of Equation (*):

    J[i+1/2] = Vdr[i+0.5]*(N[i]-(1-0.1*Pe)^5(N[i+1]-N[i])/Pe), for 0<Pe<=5;

       [i-1]  [i]  [i+1]
    --|--x--|--x--|--x--|--

    A[i]*N[i] = A[i+1]*N[i+1] + A[i-1]*N[i-1]

    Coefficients (linear geometry):
    A[i+1] = D[i+0.5]*F(|Pe[i+0.5]|)+Max(-Vdr[i+0.5];0)
    A[i-1] = D[i-0.5]*F(|Pe[i-0.5]|)+Max(Vdr[i-0.5];0)
    A[i] = A[i+1]+A[i-1]+(Vdr[i+0.5]-Vdr[i-0.5])

    where:
    F(|Pe|) = Max(0;(1-0.1*|Pe|)^5)

    In "RADIAL" case:

    AA[i+1] = A[i+1]*r[i+0.5]/(r[i]*(r[i+0.5]-r[i-0.5]))
    AA[i-1] = A[i-1]*r[i-0.5]/(r[i]*(r[i+0.5]-r[i-0.5]))
    AA[i] = AA[i+1]+AA[i-1]+(r[i+0.5]*Vdr[i+0.5]-r[i-0.5]*Vdr[i-0.5])/(r[i]*(r[i+0.5]-r[i-0.5]))

    Chemistry part considered as explicit (source term)

	//-----------------------------------------------------------------------------------------------------------------------------
	*/

	double  D_L,D_R,V_L,V_R,dL,
            Pe,fPe,
            A_L,A_R,A_C;

    double  A,B,C,F,den,
            al[LEN+1],bet[LEN+1];

    //SWEEP-SHUTTLE_CICLE*****************************************

    //Defining_sweep_coefficients
    int i;
    for(i=1;i<=LEN;i++)
    {
        //Defining_sweep_coefficients
        //Diffusion
        D_L = D_R;
        if(i==1)
            D_L = 0.5*(Di[i]+Di[i-1]);
        D_R = 0.5*(Di[i+1]+Di[i]);

        //Drift
        V_L = V_R;
        if(i==1)
            V_L = Vi[i-1];//0.5*(Vi[i]+Vi[i-1]);
        V_R = Vi[i];//0.5*(Vi[i+1]+Vi[i]);

        //DDT_Coefficients:
        //right-edge
        dL = 0.5*(l[i+2]-l[i]);
        Pe = fabs(V_R*dL/D_R);
        fPe = pow((1.0-0.1*Pe),5.0);
        //A_R = D_R*std::max(0.0,fPe)+std::max(-V_R,0.0);
        A_R = D_R*fmax(0.0,fPe)+fmax(-V_R,0.0);

        //left-edge
        dL = 0.5*(l[i+1]-l[i-1]);
        Pe = fabs(V_L*dL/D_L);
        fPe = pow((1.0-0.1*Pe),5.0);
        //A_L = D_L*std::max(0.0,fPe)+std::max(V_L,0.0);
        A_L = D_L*fmax(0.0,fPe)+fmax(V_L,0.0);

        //Accounting for Geometry Factors:
        A_R = A_R*GF_R[i];
        A_L = A_L*GF_L[i];
        A_C = A_R+A_L+GF_R[i]*V_R-GF_L[i]*V_L;

        //SWEEP Coefficients:
        A = -A_L;///[i-1](left_cell)
        B = -A_R;///[i+1](right_cell)
        C = A_C+1.0/dt;///[i](center_cell)
        F = NNi[i]*1.0/dt+Si[i];///RHS-part

        if(i==1)//see SpecTransport();
        {
            al[i-1] = al_b[0];///see_SpecTransportBoundary();
            bet[i-1] = bet_b[0];//0.0;
        }

        den = 1.0/(A*al[i-1]+C);
        al[i] = -B*den;
        bet[i] = (F-A*bet[i-1])*den;

        if(i==LEN)//see SpecTransport();
        {
            al[i] = al_b[1];///see_SpecTransportBoundary();
            bet[i] = bet_b[1];
        }
    }

    //Reverse_sweep_cycle************************************************
    for(i=LEN;i>=0;i--)
    {
        NNi[i] = al[i]*NNi[i+1]+bet[i];
        if(NNi[i]<1.e-30)
            NNi[i] = 0.0;
    }

    return &NNi[0];
}
void TransportBoundary(int N,double *Ni,double *Xi,double *Di,double *Ngas,double *Tgas,double *Tv,double *Te,double Twall)
{
    int n,nX=11;
    double Pe,D,Vt,TsrL,TsrR;
    double al,alXL,alXR,dlL,dlR,DXL,DXR,dNL,dNR;
    double Nnew,sign,NgL,NgR;

    dlL = 0.5*(l[2]-l[0]);
    dlR = 0.5*(l[LEN+2]-l[LEN]);
    //TsrL = 0.5*(Tgas[0]+Tgas[1]);
    //TsrR = 0.5*(Tgas[LEN]+Tgas[LEN+1]);
    //NgL = 0.5*(Ngas[0]+Ngas[1]);
    //NgR = 0.5*(Ngas[LEN]+Ngas[LEN+1]);

    ///Boundary_conditions********************************************

    dNL = 0.0;
    dNR = 0.0;
    for(n=0;n<N;n++)
    {
       ///Left_boundary***********************************************
        if(Gamma[n][0] == 0.0)//Ni[n][0] = Ni[n][1];
        {
            Ni[n*(LEN+2)] = Ni[n*(LEN+2)+1];
            al_bound[n][0] = 1.0;
            bet_bound[n][0] = 0.0;

            if(n==nX)
            {
                DXL = 0.5*(Di[n*(LEN+2)+1]+Di[n*(LEN+2)]);
                alXL = 1.0;
            }
        }
        else///equation for DDT with kinetic wall flux
        {
            if(Ni[n*(LEN+2)+1]-Ni[n*(LEN+2)]>0.0)
                sign = 1.0;
            else
                sign = -1.0;

            D = 0.5*(Di[n*(LEN+2)+1]+Di[n*(LEN+2)]);
            if(n==nX)
                DXL = D;

            if(n==0)
                Vt = sqrt(8*kb*Te[0]*eV_K/(pi*Mi[n]));
            else
                Vt = sqrt(8*kb*Tgas[0]/(pi*Mi[n]));

            //if(Vd!=0)
            //Vd = 0.5*(Mui[1]+Mui[0])*E;
            //Pe = Vd*0.5*(l[2]-l[0])/D;
            //al = (1.0+0.25*Gamma[n][0]*Pe*Vt/Vd)/(1.0+Pe);

            Pe = dlL*Vt/D;
            al = (1.0+sign*0.25*Gamma[n][0]*Pe);
            //al = (Tgas[0]/Tgas[1]+sign*0.25*Gamma[n][1]*Pe*TsrL/Tgas[1]);//Mankelevich_like_Jd=-D*N*dXi/dl

            if(n==nX)
                alXL = al;

            Nnew = Ni[n*(LEN+2)+1]/al;
            Ni[n*(LEN+2)] = Nnew;

            if(n!=nX)
                dNL += 0.25*Gamma[n][0]*Vt*Nnew*Mi[n];//fabs(Nnew-Nold);//

            al_bound[n][0] = 1.0/al;
            bet_bound[n][0] = 0.0;
        }

        ///Right_boundary************************************************
        if(Gamma[n][1] == 0.0)//Ni[n][LEN+1] = Ni[n][LEN];
        {
            Ni[n*(LEN+2)+LEN+1] = Ni[n*(LEN+2)+LEN];
            al_bound[n][1] = 1.0;
            bet_bound[n][1] = 0.0;

            if(n==nX)
            {
                DXR = 0.5*(Di[n*(LEN+2)+LEN+1]+Di[n*(LEN+2)+LEN]);
                alXR = 1.0;
            }
        }
        else///equation for DDT with kinetic wall flux
        {
            if(Ni[n*(LEN+2)+LEN+1]-Ni[n*(LEN+2)+LEN]>0.0)
                sign = -1.0;
            else
                sign = 1.0;

            D = 0.5*(Di[n*(LEN+2)+LEN+1]+Di[n*(LEN+2)+LEN]);
            if(n==nX)
                DXR = D;

            if(n==0)
                Vt = sqrt(8*kb*Te[LEN+1]*eV_K/(pi*Mi[n]));
            else
                Vt = sqrt(8*kb*Tgas[LEN+1]/(pi*Mi[n]));

            //if(Vd!=0)
            //Vd = 0.5*(Mui[LEN]+Mui[LEN+1])*E;
            //Pe = Vd*0.5*(l[LEN+2]-l[LEN])/D;
            //al = (1.0+0.25*Gamma[n][1]*Pe*Vt/Vd)/(1.0+Pe);//al_bound=Ni[LEN]/Ni[LEN+1]

            Pe = dlR*Vt/D;
            al = (1.0+sign*0.25*Gamma[n][1]*Pe);//al=Ni[LEN]/Ni[LEN+1]
            //al = (Tgas[LEN+1]/Tgas[LEN]+sign*0.25*Gamma[n][1]*Pe*TsrR/Tgas[LEN]);//Mankelevich_like_Jd=-D*N*dXi/dl

            if(n==nX)
                alXR = al;

            Nnew = Ni[n*(LEN+2)+LEN]/al;
            Ni[n*(LEN+2)+LEN+1] = Nnew;

            if(n!=nX)
                dNR += 0.25*Gamma[n][1]*Vt*Nnew*Mi[n];//fabs(Nnew-Nold);//

            al_bound[n][1] = al;
            bet_bound[n][1] = 0.0;
        }
    }

    ///for_N2[0]
    n = nX;
    if(dNL!=0.0)
    {
        ///Ni_version
        bet_bound[n][0] = dNL*dlL/(Mi[n]*DXL*alXL);
        Ni[n*(LEN+2)] += bet_bound[n][0];
    }

    if(dNR!=0.0)
    {
        ///Ni_version
        bet_bound[n][1] = -1.0*dNR*dlR/(Mi[n]*DXR);
        Ni[n*(LEN+2)+LEN+1] += -1.0*bet_bound[n][1]/alXR;
    }

    ///Temperature_Boundary_conditions*****************************************
    ///Left_boundary**********************
    if(!strcmp(Geom,"axial"))
        Tgas[0] = Tgas[1];
    else if(!strcmp(Geom,"cartesian"))
        Tgas[0] = Twall;

    ///Right_boundary*********************
    Tgas[LEN+1] = Twall;

    Ngas[0] = 0.0;
    Ngas[LEN+1] = 0.0;
    for(n=0;n<N;n++)
    {
        Ngas[0] += Ni[n*(LEN+2)];
        Ngas[LEN+1] += Ni[n*(LEN+2)+LEN+1];
    }

    for(n=0;n<N;n++)
    {
        Xi[n*(LEN+2)] = Ni[n*(LEN+2)]/Ngas[0];
        Xi[n*(LEN+2)+LEN+1] = Ni[n*(LEN+2)+LEN+1]/Ngas[LEN+1];
    }

    ///Ne-Te_Boundary_conditions***********************************************
    ///left_boundary**********************
    Te[0] = Te[1];
    Tv[0] = Tv[1];
    al_bound[N+1][0] = 1.0;
    bet_bound[N+1][0] = 0.0;

    ///Right_boundary**********************
    Te[LEN+1] = Te[LEN];
    Tv[LEN+1] = Tv[LEN];
    al_bound[N+1][1] = 1.0;
    bet_bound[N+1][1] = 0.0;

}
void TeTransport(double *Ne,double *Te,double *Lel,double *Vel,int Nal,double dt)
{
    double NeTe[LEN+2],De[LEN+2],Ve[LEN+1],Se[LEN+2];

    //Сетка по длине:
	/*
                left wall                                                               right wall
    Ni[i]         |                                                                       |
    Fi[i]     [0] | [1]   [2]	[3]		   		 [i-1]  [i]  [i+1]					  [I] |[I+1]
            |--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|---------->Len
    l[i]   [0]   [1]   [2]   [3]         	  [i-1]  [i]  [i+1] [i+2]              [I]  [I+1] [I+2]
    E[i]          |                                                                       |
                  |                                                                       |
    */

    NeTe[0] = 1.5*Ne[0]*Te[0];
    NeTe[LEN+1] = 1.5*Ne[LEN+1]*Te[LEN+1];

    //Defining_Drift-coefficient and Source_term******************
    int i;
    double dL;
    for(i=1;i<=LEN;i++)
    {
        NeTe[i] = 1.5*Ne[i]*Te[i];

        if(i==1)
        {
            De[0] = 2.0*Lel[0]/Ne[0]/3.0;
            De[LEN+1] = 2.0*Lel[LEN+1]/Ne[LEN+1]/3.0;
        }
        De[i] = 2.0*Lel[i]/Ne[i]/3.0;//0.666667

        dL = 0.5*(l[i+1]-l[i-1]);

        Ve[i-1] = 5.0*Vel[i-1]/3.0+0.5*(De[i]+De[i-1])*(Ne[i]-Ne[i-1])/dL;   //1.666667
        if(i==LEN)
        {
            dL = 0.5*(l[i+2]-l[i]);
            Ve[i] = 5.0*Vel[i]/3.0+0.5*(De[i]+De[i+1])*(Ne[i+1]-Ne[i])/dL;
        }

        Se[i] = 0.0;//(Jel[i]*E[i])-Ne*Sum(ki*Ni*dEei);
    }

    //проверить_геометрические факторы в поправке для сноса!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    //Transport_equation_solve************************************
    double *res;

    res = Transport_SWEEPsolve(NeTe,De,Ve,Se,&al_bound[Nal][0],&bet_bound[Nal][0],dt);
    for(i=1;i<=LEN;i++)
        Te[i] = *(res+i)/(1.5*Ne[i]);
}
<file_sep>/1D_MainFun.cpp
# include "1D_MainFun.h"

void Current_Power_calc(int,int,double *,double *,double *,double *,double *,double *,double *,double *,double *,double);
void VoltageSupply_correction(double,double,double *,double *,double *,char *,int dot);
void def_dtel_step(int,double,double *,double *,double,double,int,int,double *,double *,double *,double *);
void def_dtgas_step(int,double,double,double *,double,int,int,double *,double *,double *,double *,double *);
void mass_control(double *,double *,int,double *,double *);
double TRFaverage(double,double,double);

int main(void)
{
    extern double   Ni[Nmax][LEN+2],Xi[Nmax][LEN+2],Pgas[LEN+2],Tgas[LEN+2],
           Ngas[LEN+2],Rogas[LEN+2],Hgas[LEN+2],
           Nel[LEN+2],Te[LEN+2],Tv[LEN+2],
           Er[LEN+1],Fir[LEN+2],PowExp;

    extern double tau,dt,WextL,WextR;//dte;
    extern int N,Nneg,Npos;
    extern bool TEEDF;

    int Nedf,Nchem;
    int i,n,nt,dot,dot1,dot2,dot3,dot4,Ndot1;
    double tic,etic,period = 0.0,Wext,ticRF = 0.0;
    double te_ddt = 1.0e-12,dte = te_ddt;

    init_read();
    init_gasDBparser();
    tic = init_data();
    etic = tic;

    Mesh_gen(Len);
    Mesh_GFcalc(Geom);

    Nedf = EEDF_read_CS(N,CSFile);
    Nchem = chem_make_react(Nedf,ChemFile);
    chem_read_react(Nchem,N);

    double Vext=3000.0/Eabs,dVext=0.0,PowCalc=0.0;
    double PowTav=0.0,IcalcTav=0.0;
    double Xi_Lav[N],Ni_Lav[N],Ng_Lav,Tg_Lav,Te_Lav,Tv_Lav,Jcond_Lav;
    //double Nel_Tav[LEN+2],Te_Tav[LEN+2],Jcond_Tav[LEN+2],Er_Tav[LEN+1];
    double Kch[LEN+2][Nchem],Sch[N+1][LEN+2],Wrad[LEN+2];//Kel[LEN+2][Nedf],
    double dTe[LEN+2],dTgas[LEN+2],dNel[LEN+2],EN[LEN+2],dEzN=0.0,dEz=0.0,Icalc=0.0;
    double Del[2][LEN+2],Muel[2][LEN+2],Di[N][LEN+2],Mui[N][LEN+2],Lam[LEN+2];
    double Eold[LEN+1],Jcond[LEN+2],Jdisp[LEN+2];
    double MNT,MNTr;

    Wext = fmax(WextL,WextR);

    mass_control(&MNT,&MNTr,N,&Ni[0][0],Tgas);
    gas_LenAverage(N,Geom,&Xi[0][0],&Ni[0][0],Ngas,Tgas,Te,Tv,Xi_Lav,Ni_Lav,&Ng_Lav,&Tg_Lav,&Te_Lav,&Tv_Lav);
    gas_LenPrint(&Xi[0][0],&Ni[0][0],N,Npos,Nneg,Pgas,Tgas,Ngas,Hgas,Rogas,Nel,Te,Tv,Er,Fir,Jcond,tic);
    gas_TimePrint(Ni_Lav,N,Pgas[1],Tg_Lav,Ng_Lav,Te_Lav,Tv_Lav,0.0,Icalc,Vext,Wext,0.0,0.0,0.0,MNT,MNTr,tic,etic,dte,dt);

    for(i=1; i<=LEN; i++)
    {
        Wrad[i] = 0.0;
        Jcond[i] = 0.0;
    }

    double tg_chem,te_chem = 1e-8;
    double Ermax = 1.0e-10,Er0max = 1.0e-10;
    dot = 0, dot1 = 0, dot2 = 0, dot3 = 0, dot4 = 0, nt = 0;
    Ndot1 = 10;//int(Ndots/10);//150//20);///75);

    ///----------------------------------------------------------------------------------

    ///EEDF_Table-Data_Import/Calculation************************************************
    int err;
    err = EEDF_table_read(N,Nedf,EedfTblFl);
    if(err!=0)
        EEDF_table_calc(Nedf,N,LEN+2,&Ni[0][1],Ngas[1],Tgas[1],1.0e-5,EedfTblFl);

    ///----------------------------------------------------------------------------------

    ///*******************Main_Cycle_of_1D-RF-Discharge_Simulation*********************///
    do
    {
        dot++;
        dot1++;

        ///Ke_De_Mue_recalculation********************************************************
        for(i=1;i<=LEN;i++)
        {
            if(((nt==0) || (dot==Ndots) || (dVext>100) || (dot1==1)) && (Nel[i]>0.0))//|| || (dot1==Ndot1) || (dot1*dte/dt>0.005(dEzN>0.05)|| (dot1==1)
                EEDF_const_calc("E/N-case",Nedf,&Kch[i][0],Nel[i],Er[i],Ngas[i],&Te[i],&dTe[i],&Del[0][i],&Muel[0][i],tic);//&Kel[i][0]
                //EEDF_const_calc("Te-case",Nedf,&Kch[i][0],Nel[i],Er[i],Ngas[i],&Te[i],&dTe[i],&Del[0][i],&Muel[0][i],tic);//&Kel[i][0]

        }

        ///Transport_data_recalculation**************************************************
        if((nt==0) || (dot==Ndots)|| (dVext>100) || (dot1==1) || (dTgas[1]>10.0) || (dTe[1]*11605>100.0))///Ñðàâíåíèå òîëüêî â îäíîé òî÷êå!!!???|| (dEz/Ez>0.001)
            Trasport_coefs_calc(N,Npos+Nneg,&Ni[0][0],Er,&Del[0][0],&Muel[0][0],&Di[0][0],&Mui[0][0],Lam,Pgas,Tgas,Ngas,Te);

        ///Chem_data_Coef-s_recalculation************************************************
        for(i=1;i<=LEN;i++)
        {
            if((nt==0) || (dot==Ndots) || (dTgas[i]>30.0) || (dTe[i]*11605>100.0))
                chem_const(&Kch[i][0],Nedf,Nchem,N,Te[i],Tgas[i],tic);
        }

        ///Heat-Transport_Eq_solve*******************************************************
        HeatTransport(Hgas,Ngas,Tgas,Lam,&Ni[0][0],&Xi[0][0],&Di[0][0],N,Er,Jcond,Wrad,dt);

        ///Poisson_Equation_solve********************************************************
        //Ermax = Poisson_SORsolve(Fir,Er,Vext,&Xi[0][0],&Ni[0][0],Ngas,Npos,Nneg,tic);
        Ermax = Poisson_SORsolve(Fir,Er,Vext,&Xi[0][0],&Ni[0][0],Ngas,Npos,Nneg,etic);

        ///Chem_Rates_recalculation******************************************************
        for(i=1;i<=LEN;i++)
            chem_runge_kutta4(i,&Ni[0][i],Nneg+Npos,N,&Kch[i][0],&Sch[0][i],Nchem,&Wrad[i],dt,dte,&tg_chem,&te_chem,tic,dot);

        ///Species-DDT+Chemistry_solve***************************************************
        ///+++Charged_species:
        for(n=0;n<=Npos+Nneg;n++)
            //SpecTransport(n,&Xi[n][0],&Ni[n][0],&Di[n][0],&Mui[n][0],Ngas,Er,&Sch[n][0],dte);
            SpecTransport_expl(&Xi[n][0],&Ni[n][0],&Di[n][0],&Mui[n][0],&Sch[n][0],Ngas,Er,dte);

        ///+++Neutral_species:
        for(n=Npos+Nneg+1;n<N;n++)
            //SpecTransport(n,&Xi[n][0],&Ni[n][0],&Di[n][0],&Mui[n][0],Ngas,Er,&Sch[n][0],dt);
            SpecTransport_expl(&Xi[n][0],&Ni[n][0],&Di[n][0],&Mui[n][0],&Sch[n][0],Ngas,Er,dt);

        ///Te-equation*******************************************************************
        TeTransport_expl(Te,Nel,&Del[1][0],&Muel[1][0],Jcond,Er,&Sch[N+1][0],dTe,dte);

        ///Gas-data_recalculation********************************************************
        for(i=1;i<=LEN;i++)
            gas_TP_calc(&Ni[0][i],&Xi[0][i],N,&Pgas[i],&Tgas[i],&dTgas[i],&Ngas[i],&Rogas[i],&Hgas[i],&Nel[i],&dNel[i],&Tv[i]);

        ///Boundary-conditions_recalculation*********************************************
        TransportBoundary_mod(N,Nneg+Npos,&Ni[0][0],&Xi[0][0],&Di[0][0],&Mui[0][0],Er,Ngas,Nel,Tgas,Pgas,Tv,Te,Tw);

        ///Power_recalculation***********************************************************
        Current_Power_calc(Npos,Nneg,&Mui[0][0],&Ni[0][0],Er,Eold,Fir,Jcond,Jdisp,&Icalc,&PowCalc,dt);

        tic += dt;
        etic += dte;
        ticRF += dte;

        ///RF-Period_averaging***********************************************************
        PowTav = TRFaverage(PowCalc,Wext,dte);
        IcalcTav = TRFaverage(Icalc,Wext,dte);

        ///Writing_and_saving_Gas-data***************************************************
        if(dot==Ndots)
        {
            mass_control(&MNT,&MNTr,N,&Ni[0][0],Tgas);

            gas_LenAverage(N,Geom,&Xi[0][0],&Ni[0][0],Ngas,Tgas,Te,Tv,Xi_Lav,Ni_Lav,&Ng_Lav,&Tg_Lav,&Te_Lav,&Tv_Lav);

            gas_LenPrint(&Xi[0][0],&Ni[0][0],N,Npos,Nneg,Pgas,Tgas,Ngas,Hgas,Rogas,Nel,Te,Tv,Er,Fir,Jcond,tic);

            gas_TimePrint(Ni_Lav,N,Pgas[1],Tg_Lav,Ng_Lav,Te_Lav,Tv_Lav,Jcond_Lav,Icalc,Vext,Wext,PowCalc,PowTav,IcalcTav,MNT,MNTr,tic,etic,dte,dt);

            gas_SavePrint(&Ni[0][0],N,Pgas,Tgas,Nel,Te,Fir,tic);

            dot = 0;
        }

        ///Time-Step_Correction**********************************************************
        ///+++Charged_species:
        if((nt==0) || (dot3==1) || (fabs(1-Ermax/Er0max)>0.1))
        {
            def_dtel_step(dot3,tic,&dte,&te_ddt,te_chem,0.75,N,Npos+Nneg,Er,&Di[0][0],&Mui[0][0],&Ni[0][0]);//0.15

            Er0max = Ermax;

            dot3 = 0;
        }
        ///+++Neutral_species:
        if(nt==0 || dot2==Ndot1)
        {
            def_dtgas_step(dot2,tic,tg_chem,&dt,0.75,N,Npos+Nneg,&Xi[0][0],Rogas,Tgas,&Di[0][0],Lam);

            dot2 = 0;
        }

        ///External_Voltage-supply_correction********************************************
        if((nt==0) || (dot==Ndots) || (fabs(1.0-PowCalc/PowExp)>0.005) || (dVext>100))//&& (dot1==10)((fabs(1.0-PowTav/PowExp)>0.005) && (ticRF>1.0/Wext)))
        {
            VoltageSupply_correction(PowCalc,PowExp,&Vext,&dVext,Ngas,Geom,dot);//PowTav

            dot1 = 0;
        }

        nt++;
        dot2++;
        dot3++;

        if(ticRF>1.0/Wext)
        {
            PowTav = 0.0;
            IcalcTav = 0.0;
            ticRF = 0.0;
        }
    }
    while(tic<tau);

    return 0;
}
void Current_Power_calc(int Npos,int Nneg,double *Mui,double *Ni,double *E,double *Eold,double *Fi,double *Jcond,double *Jdisp,double *Icalc,double *PowCalc,double dt)
{
    int i,n;
    double Pow[LEN+2];

    for(i=1;i<=LEN;i++)
    {
        ///Conductivity_Current****************************
        Jcond[i] = Ni[i]*Mui[i];///electrons
        for(n=1; n<=Npos; n++) ///positive_ions
            Jcond[i] += Ni[n*(LEN+2)+i]*Mui[n*(LEN+2)+i];
        for(n=Npos+1; n<=Npos+Nneg; n++) ///negative_ions
            Jcond[i] += Ni[n*(LEN+2)+i]*Mui[n*(LEN+2)+i];
        Jcond[i] *= e*E[i];//[Abs/cm2*s]
        Jcond[i] = fabs(Jcond[i]);//[Abs/cm2*s]

        ///Displacement_Current****************************
        Jdisp[i] = (E[i]-Eold[i])/(dt*4.0*pi)*Eabs;

        ///Power_distribution******************************
        //Pow[i] = Jcond[i]*fabs(Fi[i]);
        Pow[i] = Jcond[i]*fabs(E[i]);///[abs/cm^3]
        //Pow[i] = (Jcond[i]+fabs(Jdis[i]))*fabs(Fi[i]);
    }

    ///Summarized_Current**********************************
    *Icalc = eKl/e*LENintegral(Jcond,Geom);//[Kl/s=A]

    ///Power_Supply_calculation****************************
    *PowCalc = eKl/e*LENintegral(Pow,Geom)*Eabs/Hght;///[W/cm^2]
}
void VoltageSupply_correction(double PowCalc,double PowExp,double *Vext,double *dVext,double *Ngas,char *Geom,int dot)
{
    int i;
    double dV,Vnew,
        Vmax=10000/Eabs;

    if(*dVext>100.0)
        *dVext = 0.0;

    double th = 1.05;
    double Vmin = 10.0;///[V]

    Vnew = *Vext;///[abs]
    if(fabs(1.0-(PowCalc/PowExp))>0.005)
    {
        if(PowCalc>PowExp)
            Vnew /= 1.005;
        else
            Vnew *= 1.005;
    }
    if(PowCalc/PowExp>th)
        Vnew /= th;

    ///Restriction_on_Vext*******************************************
    if(Vnew>Vmax)///[abs]
        Vnew = Vmax;

    if((Vnew!=*Vext) && dot==Ndots)
        printf("Vext = %.1lf[V]\n",Vnew);

    *dVext = fabs(Vnew - *Vext);
    *Vext = Vnew;

    /*
    double res;
    dEN = 0.0;
    if(!strcmp(Geom,"axial"))
    {
        for(i=1; i<=LEN; i++)
        {
            res = fabs(EN[i] - Enew/Ngas[i]*1e17*Eabs);//[Td]
            EN[i] = Enew/Ngas[i]*1e17*Eabs;//[Td]

            dEN += dEN*0.5*(l[i+1]+l[i])*(l[i+1]-l[i]);
        }
        dEN *= 2.0/(Len*Len);
    }
    else if(!strcmp(Geom,"cartesian"))
    {
        for(i=1; i<=LEN; i++)
        {
            res = fabs(EN[i] - Enew/Ngas[i]*1e17*Eabs);//[Td]
            EN[i] = Enew/Ngas[i]*1e17*Eabs;//[Td]

            dEN += dEN*(l[i+1]-l[i]);
        }
        dEN /= Len;
    }

    *dEzN += dEN;
    */
}
void def_dtel_step(int dot,double tic,double *dte,double *te_min0,double te_chem,double CFL,int N,int Nion,double *E,double *Di,double *Mui,double *Ni)
{
    double te_min, te_drift, te_diff=1.0, te_ddr=1.0;
    double dl,Cmax;

    Cmax = 0.15;//0.75;////0.75

    ///electronic_time_step***********************************************************
    te_drift = (l[1]-l[0])/fabs(0.5*(Mui[0]+Mui[1])*E[0]);//only_electrons
    for(int i=1; i<=LEN; i++)
    {
        dl = l[i+1]-l[i];

        te_diff = fmin(dl*dl/(2*Di[i]),te_diff);//only_electrons

        te_drift = fmin(dl/(fabs(0.5*(Mui[i+1]+Mui[i])*E[i])),te_drift);//only_electrons
        //te_drift = fmin(dl/(fabs(0.5*(Mui[i+1]+Mui[i])*E[i]*Ni[i]/Ni[0])),te_drift);//weighted

        te_ddr = fmin(fabs(dl/(2*Di[i]/dl-(0.5*(Mui[i+1]+Mui[i])*E[i]))),te_ddr);
    }

    te_min = fmin(te_diff,te_drift);
    te_min = fmin(te_ddr,te_min);
    te_min *= Cmax*CFL;
    //if(te_chem>=1.e-10)
        te_min = fmin(te_min,te_chem);

    ///dte_time-step_correction*******************************************************
    if(*dte<te_min)// && (Ermax/Er0max<0.99))//(fabs(dEr-dEr0)<0.1));//(dEr/dEr0<0.999))////(dEr/dEr0<0.95))//(fabs(dEr-dEr0)<0.1))//))(dEr<dEr0))
        *dte *= 1.01;
    if(te_min/(*te_min0)>=1.0001)// && (Ermax/Er0max<0.99))//(fabs(dEr-dEr0)<0.1));//(dEr/dEr0<0.999))////(dEr/dEr0<0.95))//(fabs(dEr-dEr0)<0.1))//))(dEr<dEr0))
        *dte *= 1.001;
    if(te_min/(*te_min0)<0.999 && *dte>te_min)
        *dte /= 1.005;//*te_min0/te_min;
    //if((1.e-12<dte<te_min) && (Ermax/Er0max>=1.01))//dEr/dEr0>=1.005))//(dEr-dEr0)>=0.05))//(dEr/dEr0>1.01))//(fabs(dEr-dEr0)>=0.1))//(dEr>dEr0))
    //dte *= 0.99;//1.005*dEr/dEr0;//*dEr/dEr0;//0.9975;
    //if((1.e-12<dte<te_min) && (Ermax/Er0max>=1.05))//dEr/dEr0>=1.005))//(dEr-dEr0)>=0.05))//(dEr/dEr0>1.01))//(fabs(dEr-dEr0)>=0.1))//(dEr>dEr0))
    //dte *= 0.99*Er0max/Ermax;//1.005*dEr/dEr0;//*dEr/dEr0;//0.9975;
    //if(dte>=0.1*te_min)
    //dte *= 0.999;
    if(*dte<1.e-12)
        *dte = 1.e-12;
    if(*dte>1.e-9)
        *dte = 1.e-9;
    if(*dte>te_min)
        *dte=te_min;

    *te_min0 = te_min;

    //*dte = 5.e-12;

    ///Logging************************************************************************
    FILE *Tmin;

    if(dot==0)
    {
        Tmin = fopen("Dte.txt", "w");
        fprintf(Tmin,"t,[s]\tdte\tte_min\tte_diff\tte_drift\tte_chem\tte_ddr\n");
    }
    else
        Tmin = fopen("Dte.txt", "a+");

    fprintf(Tmin,"%.10e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",tic,*dte,*te_min0,te_diff,te_drift,te_chem,te_ddr);
    fclose(Tmin);
}
void def_dtgas_step(int dot,double tic,double tg_chem,double *dtg,double CFL,int N,int Nion,double *Xi,double *Rogas,double *Tgas,double *Di,double *Lam)
{
    double tg_diff=1.0, tg_therm=1.0, tg_min;
    double Dmax,Cp,dl,Th;
    int i,ni,n;

    ///gaseous_time_step*************************************************************
    for(i=1; i<=LEN; i++)
    {
        dl = l[i+1]-l[i];
        gas_HCpSi_calc(Tgas[i],N);

        Cp = 0.0;
        Dmax = 0.0;
        for(n=1; n<N; n++)
        {
            ni = n*(LEN+2)+i;

            if(n>Nion)
                Dmax = fmax(Dmax,Di[ni]);

            Cp += HCpSi[1][n]*Xi[ni];//[ýðã/ã*Ê]
        }
        ///diffusional
        tg_diff = fmin(dl*dl/(2*Dmax),tg_diff);

        ///thermal
        Th = Lam[i]/(Rogas[i]*Cp);//Lam-[Ýðã/(ñ*ñì*K)]
        tg_therm = fmin(dl*dl/(2*Th),tg_therm);

    }
    tg_min = CFL*fmin(tg_therm, tg_diff);
    tg_min = fmin(tg_min,tg_chem);//tg_chem
    //tg_min *= CFL;

    ///dt_time-step_correction*******************************************************
    if(*dtg<tg_min)// && (Ermax/Er0max<0.99))//(fabs(dEr-dEr0)<0.1));//(dEr/dEr0<0.999))////(dEr/dEr0<0.95))//(fabs(dEr-dEr0)<0.1))//))(dEr<dEr0))
        *dtg *= 1.01;
    if(tg_min/(*dtg)<0.95)
        *dtg /= 1.05;//*te_min0/te_min;
    if(*dtg>=tg_min)
        *dtg=0.99*tg_min;
    if(*dtg<1.e-10)
        *dtg = 1.e-10;

    ///Logging************************************************************************
    FILE *Tmin;

    if(dot==0)
    {
        Tmin = fopen("Dt_gas.txt", "w");
        fprintf(Tmin,"t,[s]\tdt_gas\ttg_min\ttg_diff\ttg_therm\ttg_chem\n");
    }
    else
        Tmin = fopen("Dt_gas.txt", "a+");

    fprintf(Tmin,"%.7e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",tic,*dtg,tg_min,tg_diff,tg_therm,tg_chem);

    fclose(Tmin);
}
void mass_control(double *MNT,double *MNTr,int N,double *Ni,double *Tgas)
{
    int i,n,ni;
    double Sum,Sumr;

    Sum = 0.0;
    Sumr = 0.0;
    for(i=0; i<=LEN+1; i++)
    {
        for(n=1; n<N; n++)
        {
            ni = n*(LEN+2)+i;

            Sum += Mi[n]*Ni[ni];//*Tgas[i];
            Sumr += Mi[n]*Ni[ni]*(i+0.5);//*Tgas[i];
        }
    }
    Sumr /= (LEN+1)*(LEN+1);

    *MNT = Sum;
    *MNTr = Sumr;
}
bool converge(double *Y0, double *Y, double eps)
{
    double Norm, cnv;
    bool Conv;
    int i;

    cnv = 0.0;
    for(i=1; i<=LEN; i++)
    {
        //Norm = (Y[i]-Y0[i])*(Y[i]-Y0[i]);
        //Norm = sqrt(Norm);

        Norm = fabs(1-Y[i]/Y0[i]);

        if(Norm<eps)//(fabs(Fi[i]-Res)<exact);//1.0-Res/Fi[i])<1.0e-5
            cnv++;

        Y0[i] = Y[i];
    }

    if(cnv<LEN)
        Conv = false;
    else
        Conv = true;

    return Conv;
}
int sign(double x)
{
    int y;

    if(x>0)
        y = 1;
    else if(x<0)
        y = -1;
    else
        y = 0;

    return y;
}
double LENintegral(double *Y,char *Geom)
{
    int i;
    double Yy = 0.0;

    if(!strcmp(Geom,"axial"))
    {
        for(i=1; i<=LEN; i++)
            Yy += Y[i]*0.5*(l[i+1]+l[i])*(l[i+1]-l[i]);
        Yy *= 2*pi;//[ÑÃÑ/ñ]
    }
    else if(!strcmp(Geom,"cartesian"))
    {
        for(int i=1; i<=LEN; i++)
            Yy += Y[i]*(l[i+1]-l[i]);
        Yy *= Hght;//[ÑÃÑ/ñ]
    }

    return Yy;
}
double LENaverage(double *Y,char *Geom)
{
    double Yy;

    if(!strcmp(Geom,"axial"))
        Yy = LENintegral(Y,Geom)/(pi*Len*Len);
    else if(!strcmp(Geom,"cartesian"))
        Yy = LENintegral(Y,Geom)/(Len*Hght);

    return Yy;
}
double TRFaverage(double Y,double Wext,double dt)
{
    double Yav;

    Yav += Y*dt*Wext;

    return Yav;
}
<file_sep>/1DPoisson.h
//Защита от повторного включения заголовка
#ifndef POISSON_CALC_H
#define POISSON_CALC_H

double Poisson_SORsolve(double *,double *,double,double *,double *,double *,int,int,double);
void Poisson_boundary(int,double *,int,double,double,double,double,double);

#endif
<file_sep>/Initials.cpp
# include "1D_MainFun.h"
# include "Initials.h"

double Xi[Nmax][LEN+2],Ni[Nmax][LEN+2],Mi[Nmax],LJi[Nmax][2],Pgas[LEN+2],Tgas[LEN+2],Ngas[LEN+2],Rogas[LEN+2],Hgas[LEN+2];
double Nel[LEN+2],Te[LEN+2],Tv[LEN+2];
double Gamma[Nmax][2];
double Ez,Er[LEN+1],Fir[LEN+2],PowExp;
double Tinit,Pinit,Xinit[Nmax],E_Ninit;
double dTgas,dTe,dNel;
double Len,Hght,Tw;
double tau,dt;
int vL,vR;
double VolL,VolR,EpsL,EpsR,DepL,DepR,WextL,WextR;
int N,NR,Nt,Nte,Ndots,Nneg,Npos;
char Spec[Nmax][10],Spec_R[Nmax][10],GeomVect[10];
char CSFile[50],ChemFile[50],EedfTblFl[50];
char Geom[10],Start[10],FileStart[50];
double CXi[Nmax][2][8];
double CDi[Nmax][25],CMui[Nmax];
bool TEEDF;

void init_read()//считывание начальных данных
{
    char symb[] = "------------------------------------------------------------";

	FILE *init;
	//init = fopen("Comps_N2.txt", "r");
	//init = fopen("Comps_Ar.txt", "r");
    init = fopen("Comps_O2.txt", "r");

    char Cmt[300];
    //считывание комментария
	fscanf(init,"%s",&Cmt);
	do
	{
		fscanf(init,"%s",&Cmt);
	}while(strcmp(Cmt,symb)!=0);

	//считывание начальных условий, названий компонент и их концентрации*************************
	///*******************************Initial_data***********************
	fscanf(init,"%s",&Cmt);
	fscanf(init,"%d%s",&N,&Cmt);
	fscanf(init,"%lf%s",&tau,&Cmt);
	fscanf(init,"%lf%s",&dt,&Cmt);
	fscanf(init,"%d%s",&Ndots,&Cmt);
	fscanf(init,"%lf%s",&Pinit,&Cmt);
	fscanf(init,"%lf%s",&Tinit,&Cmt);
	fscanf(init,"%lf%s",&Len,&Cmt);
	fscanf(init,"%lf%s",&Hght,&Cmt);
	fscanf(init,"%lf%s",&Tw,&Cmt);
	fscanf(init,"%lf%s",&PowExp,&Cmt);

    ///*****************************Additional_data**********************
	fscanf(init,"%s",&Cmt);
	fscanf(init,"%s%s",&Geom,&Cmt);
	fscanf(init,"%s%s%s",&Start,&FileStart,&Cmt);
	fscanf(init,"%s%s",&CSFile,&Cmt);
	fscanf(init,"%s%s",&ChemFile,&Cmt);

	///Boltzmann-calc_or_Table
	char EedfCase[10];
        fscanf(init,"%s%s",&EedfCase,&Cmt);
    if(!strcmp(EedfCase,"stable"))
        TEEDF = true;
    else if(!strcmp(EedfCase,"recalc"))
        TEEDF = false;
    fscanf(init,"%s%s",&EedfTblFl,&Cmt);

    ///*************************Voltage-Supply***************************
    fscanf(init,"%s",&Cmt);
    fscanf(init,"%s",&Cmt);
	fscanf(init,"%d%lf%lf%lf%lf",&vL,&VolL,&WextL,&EpsL,&DepL);
	fscanf(init,"%d%lf%lf%lf%lf",&vR,&VolR,&WextR,&EpsR,&DepR);
	VolL /= Eabs;
	VolR /= Eabs;
	WextL *= 1.e6;///[MHz]
	WextR *= 1.e6;///[MHz]

    ///*******************************Species_data***********************
	fscanf(init,"%s",&Cmt);
	int n = 0,N1;
	Npos = 0;
	Nneg = 0;
	fscanf(init,"%s",&Cmt);///Form
	fscanf(init,"%d%s%lf%lf%lf",&N1,&Spec[n],&Xinit[n],&Gamma[n][0],&Gamma[n][1]);///electrons
	for(n=1;n<N;n++)
    {
        fscanf(init,"%d%s%lf%lf%lf",&N1,&Spec[n],&Xinit[n],&Gamma[n][0],&Gamma[n][1]);
        if(strchr(Spec[n],'+')!=NULL)
            Npos++;
        else if(strchr(Spec[n],'-')!=NULL)
            Nneg++;
    }

  	///**********************Print_Rates_for_components******************
	fscanf(init,"%s",&Cmt);
	NR = 0;
	int err;
	fscanf(init,"%s",&Spec_R[NR]);
	while(strcmp(Spec_R[NR],symb)!=0)
	{
		err = 1;
		for(n=0;n<N;n++)
		{
			if(!strcmp(Spec[n],Spec_R[NR]))
			{
				NR += 1;
				err = 0;
				break;
			}
		}

		if(err == 1)
			printf("!!Warning:Unknown component - %s for ChemRates analysis!!Excluded!\n\n",Spec_R[NR]);

		fscanf(init,"%s",&Spec_R[NR]);
	}
	fclose(init);
}
double init_data()//задание начальных условий
{
	int i,n,k,st;
	double tic;

	if(!strcmp(Start,"saved"))
    {
        char Cmt[20];
        FILE *save;

        save = fopen(FileStart, "r");

        //Парсинг_файла***************************************************

        fscanf(save,"%s%lf",&Cmt,&tic);///Time
        fscanf(save,"%s%lf",&Cmt,&Pinit);///Pressure_[Torr]
        fscanf(save,"%s%lf",&Cmt,&Ez);///Field_[V/cm]

        fscanf(save,"%s",&Cmt);///Temperature_[K]
        for(i=0;i<=LEN+1;i++)
            fscanf(save,"%lf",&Tgas[i]);

        fscanf(save,"%s",&Cmt);///Electronic_Temperature_[eV]
        for(i=0;i<=LEN+1;i++)
            fscanf(save,"%lf",&Te[i]);

        fscanf(save,"%s",&Cmt);///Field(radial)_[СГС]
        for(i=0;i<=LEN+1;i++)
        {
            fscanf(save,"%lf",&Fir[i]);//[Volt]
            Fir[i] /= Eabs;
        }

        for(n=0;n<N;n++)///Ni_concentration_[cm-3]
        {
            fscanf(save,"%s",&Cmt);
            for(i=0;i<=LEN+1;i++)
                fscanf(save,"%lf",&Ni[n][i]);
        }

        fclose(save);

        st = 1;
    }
    else
        st = 0;

	//Газовые компоненты**********************************************
	for(i=0;i<=LEN+1;i++)
    {
        Pgas[i] = Pinit*p0; //Torr

        if(st==0)
            Tgas[i] = Tinit-i*i*(Tinit-Tw)/(LEN+1)/(LEN+1);

        gas_HCpSi_calc(Tgas[i],N);

        Ngas[i] = Pgas[i]/(kb*Tgas[i]);
        Hgas[i] = 0.0;
        Rogas[i] = 0.0;
        for(n=0;n<N;n++)
        {
            if(st==0)
            {
                Xi[n][i] = Xinit[n];
                Ni[n][i] = Xinit[n]*Ngas[i];
            }
            else
                Xi[n][i] = Ni[n][i]/Ngas[i];

            Rogas[i] += Ni[n][i]*Mi[n];
            if(n>=1)
                Hgas[i] += HCpSi[0][n]*Ni[n][i]*Mi[n];//[эрг/cm^3]
        }
        Nel[i] = Ni[0][i];
    }
	//****************************************************************

	//Потенциал электрического поля***********************************

	if(st==0)
        Ez = E_Ninit*Ngas[0]*1e-17;//E_N[Td] = E[B/cm]*1e17/Ngas[cm-3];
    Ez /= Eabs;//E=E[B/cm]/E[abs]

    if(st==0)
    {
        for(i=0;i<=LEN;i++)
        {
            Fir[i] = 0.0;
            Er[i] = 0.0/Eabs;
        }
        Fir[LEN+1] = 0.0;
    }

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

	//Начальная Te****************************************************
    for(i=0;i<=LEN+1;i++)
    {
        if(st==0)
            Te[i] = 0.5;
    }
	//****************************************************************

	//Writing_data****************************************************
	init_print();
	//****************************************************************

	if(st==0)
        tic = 0.0;

	return tic;
}
void init_print()//запись начальных данных
{
	char symb[] = "------------------------------------------------------------";

	printf("%s\nInitial Data:\n\n",symb);
	printf("t = %.2lf[s]\n",0.0);
	printf("I = %.2lf[W]\tXe = %.2e\tTe = %.1f[eV]\n",PowExp,Nel[1]/Ngas[1],Te[1]);
	printf("P = %.1f[Torr]\tT = %.1f[K]\t%s = %.2e\n",Pgas[1]/p0,Tgas[1],Spec[0],Ni[0][1]);
	printf("\n%s\n\n",symb);

	//File_rewriting**************************************************

    FILE *log;
    log = fopen("Gas_data.txt", "w");
	fclose(log);
}
void init_gasDBparser()//получение параметров газовых компонент из базы данных
{
    char symb[] = "------------------------------------------------------------";
    char symb2[] = "//----------------";

    char Cmt[300];
    char str[10],strv[10],VCmt[100];
    double we,xewe,dH0,Enull;

    FILE *db;
	db = fopen("GasPropDB.txt", "r");

    int n,x,i,k,Nerr,err;
    int nvib = 0;
    err = 0;
    Nerr = 0;
    Mi[0] = me;//electrons
    for(n=1;n<N;n++)
    {
        //проверка элемента:колебательно-возбужденное состояние или нет
        err = 0;
        i = 0;
        while((Spec[n][i]!='(') && (Spec[n][i]!='\0'))
            i++;

        if(Spec[n][i]=='(')
        {
            memset(str, 0, sizeof(str));
            int k = 0;
            while(Spec[n][k]!='(')
            {
                str[k] = Spec[n][k];
                k++;
            }
            nvib = 1;

            memset(VCmt, 0, sizeof(VCmt));
            sprintf(VCmt,"//%s_Spectrum_data:",str);
        }
        else
        {
            memset(str, 0, sizeof(str));
            strcat(str,Spec[n]);
        }

        //поиск_по_базе_данных**************************
        while(strcmp(Cmt,str)!=0)
        {
            fscanf(db,"%s",&Cmt);

            if(!strcmp(Cmt,"END"))
            {
                err++;
                break;
            }
        }

        if(err==0)
        {
            ///Gas_Parameters_Reading****************************************************
            fscanf(db,"%lf%lf%lf",&Mi[n],&LJi[n][0],&LJi[n][1]);
            Mi[n] *= ma;
            for(x=0;x<2;x++)
            {
                for(i=0;i<8;i++)
                    fscanf(db,"%lf",&CXi[n][x][i]);//CXi[n][x][7]-доп.энтальпия образования компонент от основного состояния (не учт. в разложении)
            }

            ///Transport_Data_Reading****************************************************
            fscanf(db,"%s%lf%lf",&Cmt,&CDi[n][0],&CDi[n][1]);
            if((strchr(Spec[n],'+')!=NULL) || (strchr(Spec[n],'-')!=NULL))///Ions
            {
                if(strchr(Spec[n],'+')!=NULL)
                    CMui[n] = 1;
                else if(strchr(Spec[n],'-')!=NULL)
                    CMui[n] = -1;

                CDi[n][2] = 0.0;///Number_of_points
                if((strchr(Cmt,'E/N')!=NULL))
                {
                    fscanf(db,"%s",&Cmt);
                    k=3;
                    fscanf(db,"%s",&Cmt);
                    do
                    {
                        CDi[n][k] = atof(Cmt);
                        fscanf(db,"%lf",&CDi[n][k+1]);
                        k+=2;
                        fscanf(db,"%s",&Cmt);
                    }while(strcmp(Cmt,symb2)!=0);
                    CDi[n][2] = int((k-3)/2);///Number_of_points
                }
            }

            ///Drift-Diff_Coefs_correction**************
            double  DPT,MuPTi;
            MuPTi = 1.0;///Mankelevich???
            DPT = 760*p0/pow(273.15,CDi[n][1]);
            if(n<Npos+Nneg)///ions
            {
                DPT *= pow(273.15/300.0,CDi[n][1]);
                CDi[n][0] *= DPT;
                CMui[n] *= MuPTi;
            }
            else///neutrals
                CDi[n][0] *= DPT;

                //"N2+,N4+" - CDi=0.07,0.058;[cm2/s]--из справочника "Физ.Величины"_стр.433
                //"[N,N(2D),N(2P),N2..." - CDn = 0.06,0.028;[cm2/s]--MankModel-data

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


            ///Vibrational_Spectrum_Data*************************************************
            if(nvib>0)
            {
                while(strcmp(Cmt,VCmt)!=0)
                    fscanf(db,"%s",&Cmt);
                fscanf(db,"%lf%lf",&we,&xewe);

                Enull = (0.5*we-0.5*0.5*xewe)/cm_eV;

                nvib -= 1;
                memset(strv, 0, sizeof(strv));
                sprintf(strv,"%s(%d)",str,nvib);
                while(!strcmp(strv,Spec[n+nvib]))
                {
                    dH0 = ((nvib+0.5)*we-(nvib+0.5)*(nvib+0.5)*xewe)/cm_eV-Enull;

                    Mi[n+nvib] = Mi[n];
                    LJi[n+nvib][0] = LJi[n][0];
                    LJi[n+nvib][1] = LJi[n][1];
                    for(x=0;x<2;x++)
                    {
                        for(i=0;i<8;i++)
                            CXi[n+nvib][x][i] = CXi[n][x][i];
                        if(nvib>0)
                            CXi[n+nvib][x][7] += dH0;
                    }
                    CDi[n+nvib][0] = CDi[n][0];
                    CDi[n+nvib][1] = CDi[n][1];

                    nvib ++;

                    sprintf(strv,"%s(%d)",str,nvib);
                }

                n += nvib-1;
                nvib = 0;
            }

            rewind(db);
        }
        else
        {
            printf("ERROR!!!_Element %s was not found in the DataBase!!!\n", Spec[n]);
            Nerr ++;

            rewind(db);
        }

    }

    fclose(db);

    //Запись_лога**************************************************************
    FILE *log;
    log = fopen("Log_GasDB.txt", "w");

    fprintf(log,"Parsed elements' parameters from LogGasDB.txt\n\n");
    for(n=1;n<N;n++)
        fprintf(log,"%d\t%s\t%.3lf\t%.2lf\t%.2lf\n",n,Spec[n],Mi[n]/ma,LJi[n][0],CXi[n][0][7]);

    fclose(log);
    //*************************************************************************

    //выдача отчета об ошибках:
    if(Nerr==0)
        printf("\nAll Components were parsed from DataBase\n");
    else
    {
        char a;
        printf("\nPlease, fill the DataBase!!Nerr=%d\n",Nerr);
        scanf("%c",&a);
        exit(1);
    }
}
<file_sep>/EEDF_calc.cpp
# include "1D_MainFun.h"
# include "EEDF_calc.h"

# define NEmax 3000

int CStype = 10,CSout;
int CSref[Nmax][10][CSmax],CSR[CSmax][5];
double CS[NEmax][CSmax],Ith[CSmax];
double EedfTbl[300][CSmax+7];
double EN_Tbl,dEN_Tbl;
double Emax = 50.0;
double dEev = Emax/NEmax,
       dE = dEev*1.602e-12;//[eV]=1.602e-12[erg]

int EEDF_read_CS(int N,char *CSFile)//считывание сечений EEDF-процессов(возврат кол-ва реакций)
{
	FILE *cross;
    cross = fopen(CSFile, "r");

	FILE *log;
	log = fopen("Log_CS.txt", "w");
	fclose(log);

	log = fopen("Log_CS.txt", "a+");

	double X[150],Y[150],A[150],B[150];
	char CSstr[100],CSstr1[10],Cmt[100];
	double Stat,shift;

	char KeyW[][20] =
	{
		"ELASTIC",
		"ROTATION",
		"VIB-EXCITATION",
		"EXCITATION",
		"DISSOCIATION",
		"DEEXCITATION",
		"IONIZATION",
		"ATTACHMENT",
		"RECOMBINATION",
		"DISSATTACHMENT"
	};

	char symb[] = "------------------------------------------------------------";

	int i,k,K,j,J,err,n,n_t,Kw,add,d,s,Serr,rev;

	//считывание комментария
	fscanf(cross,"%s",&Cmt);
	do
	{
		fscanf(cross,"%s",&Cmt);
	}while(strcmp(Cmt,symb)!=0);

	j = 0;
	Serr = 0;
	fscanf(cross,"%s",&CSstr);//считывание первой строки
	while(strcmp(CSstr,"END")!=0)
	{
		add = 0;
		rev = 0;

		err = 0;
		for(i=0;i<CStype;i++)//сравнение с ключевыми словами
		{
			if(!strcmp(CSstr,KeyW[i]))
			{
				Kw = i;
				CSR[j][0] = Kw;//CS-type

				err = 1;
				break;
			}
		}
		if(err==0)//нет совпадений
		{
			printf("!!!Can't find Keyword!!! %s\n", CSstr);
			fprintf(log,"!!!Can't find Keyword!!! %s\n", CSstr);

			Serr += 1;
			break;
		}

		//считывание второй строки
		fscanf(cross,"%s",&CSstr);
		err = 0;
		for(n=0;n<N;n++)
		{
			if(!strcmp(CSstr,Spec[n]))
			{
				n_t = n;///target particle
				CSR[j][2] = n_t+1;
				err += 1;
				break;
			}
		}
		if(err==0)
		{
			printf("!!!Unknown target particle - %s in #%d cross section!!!\n ",CSstr,j);
			fprintf(log,"!!!Unknown target particle - %s in #%d cross section!!!\n ",CSstr,j);

			Serr += 1;
			///break;
		}

		///обработка 2-й строки
		if(Kw>=1)//w/o ELASTIC+ROT-EXCITATION
		{
			fscanf(cross,"%s%s",&Cmt,&CSstr);
			if(Kw==4 || Kw==9)
				fscanf(cross,"%s",&CSstr1);

			if(!strcmp(Cmt,"="))//forvard or rev
            {
                CSR[j][1] = -1;
				rev = 1;

				fscanf(cross,"%lf",&Stat);///Stat=g_gr/g_ex=st_w[nr_p]/st_w[nr_t];
            }
			if(!strcmp(Cmt,"->"))
                CSR[j][1] = 1;

			err = 0;
			for(n=0;n<N;n++)
			{
				if(!strcmp(CSstr,Spec[n]))
				{
					CSR[j][3] = n+1;//product particle
					err += 1;

					if(err==2)
						break;
				}
				if((!strcmp(CSstr1,Spec[n]))&&(Kw==4 || Kw==9))///"DISSOCIATION" OR "DISSATTACHMENT"
				{
					CSR[j][4] = n+1;//product particle
					err += 1;

					if(err==2)
						break;
				}
			}

			///предупреждение
			if((Kw==4 || Kw==9)&&(err<2))
			{
				printf("!!!Unknown products in CS: %s %s %s %s - !!!\n ",Spec[n_t],Cmt,CSstr,CSstr1);
				fprintf(log,"!!!Unknown products in CS: %s %s %s %s - !!!\n ",Spec[n_t],Cmt,CSstr,CSstr1);

				Serr += 1;
				///break;
			}
			else if((err==0)&&(Kw!=4 || Kw!=9))
			{
				printf("!!!Unknown product particle in CS: %s %s %s - !!!\n ",Spec[n_t],Cmt,CSstr);
				fprintf(log,"!!!Unknown product particle in CS: %s %s %s - !!!\n ",Spec[n_t],Cmt,CSstr);

				Serr += 1;
				///break;
			}
		}

		///обработка дополнения ко 2-й строке и считывание 3-ей строки
		fscanf(cross,"%s",CSstr);
		if(!strcmp(CSstr,"//"))//доп. реакции с таким же сечением "//add: N2(3) N2(5) //"
		{
			add = 0;
            fscanf(cross,"%s",CSstr);
			while(strcmp(CSstr,"//")!=0)
			{
				err = 0;
				for(n=0;n<N;n++)
				{
					if(!strcmp(CSstr,Spec[n]))
					{
						add += 1;
						CSR[j+add][0] = CSR[j][0];//type
						CSR[j+add][1] = CSR[j][1];//forvard or rev
						CSR[j+add][2] = n+1;//target particle
						CSR[j+add][3] = CSR[j][3];//product particle
						if(Kw==4)
                            CSR[j+add][4] = CSR[j][4];//product particle
						err += 1;
						break;
					}
				}

				if(err==0)
				{
					printf("!!!Unknown product particle in CS: %s -> add: %s - !!!\n ",Spec[n_t],CSstr);
					fprintf(log,"!!!Unknown product particle in CS: %s -> add: %s - !!!\n ",Spec[n_t],CSstr);

					Serr += 1;
					break;
				}

				fscanf(cross,"%s",CSstr);
			}

            ///сдвиг сечения или нет
            shift = 0;
			fscanf(cross,"%s",CSstr);
			if(!strcmp(CSstr,"shift"))
                shift = 1;

			fscanf(cross,"%lf%s",&Ith[j],Cmt);

			for(d=1;d<=add;d++)
			{
				int n0;
				n0 = CSR[j][2]-1;
				n = CSR[j+d][2]-1;
				if(shift>0)
                    Ith[j+d] = Ith[j] - (CXi[n][0][7]-CXi[n0][0][7]);///same product case!!
                else
                    Ith[j+d] = Ith[j];///без сдвига
			}
		}
		else
		{
			Ith[j] = atof(CSstr);
			if(Kw==1)//rotational
				Ith[j] = Ith[j]/11605;//[eV]

			fscanf(cross,"%s",Cmt);
		}

		///формирование массива ссылок на номер сечения
		for(d=0;d<=add;d++)
		{
			n_t = CSR[j+d][2]-1;

			s = CSref[n_t][Kw][0]+1;
			CSref[n_t][Kw][s] = j+d+1;//j+1??????????????????****************
			CSref[n_t][Kw][0] += 1;
		}

		///считывание 4-ей строки
		fscanf(cross,"%s",Cmt);

		///считывание и обработка сечения
		fscanf(cross,"%s",CSstr);
		if(!strcmp(CSstr,symb))
		{
			k = 0;
			//X[0] = 0.0;
			//Y[0] = 0.0;
			fscanf(cross,"%s",CSstr);
			while(strcmp(CSstr,symb)!=0)
			{
				X[k] = atof(CSstr);
				fscanf(cross,"%lf",&Y[k]);
				k += 1;
				fscanf(cross,"%s",CSstr);
			}
			K = k; //счётчик
			X[K] = 1500;//Emax;
			Y[K] = 0.0;//

			///Loging*****************************************************
			fprintf(log,"#CS=%d %s for particle %s\n",j,KeyW[Kw],Spec[n_t]);
			for(k=0;k<K;k++)
				fprintf(log,"%.2lf\t%.2e\n",X[k],Y[k]);
			fprintf(log,"\n");
			///Loging*****************************************************
		}
		else
		{
			printf("WARNING!Error while reading CS %s\t for particle %s\n",KeyW[Kw],Spec[n_t]);
			fprintf(log,"WARNING!Error while reading CS %s\t for particle %s\n",KeyW[Kw],Spec[n_t]);
			Serr += 1;
			break;
		}

		///линеаризация сечения
		for(k=0; k<K; k++)
		{
			A[k] = (Y[k+1]-Y[k])/(X[k+1]-X[k]);
			B[k] = Y[k]-A[k]*X[k];
		}
		//экстраполяция до лимита по энергии - сохранение наклона!!!!!!!!
		//A[K-1] = (A[K-2]+A[K-3]+A[K-4])/3;
		//B[K-1] = B[K-2];//(B[K-2]+B[K-3]+B[K-4])/3;

		///построение сечений по расчётной энергетической сетке
		k = 0;
		double Ei, CSj;
		for(i=0; i<NEmax; i++)
		{
			Ei = (i+0.5)*dEev;
			if(Ei<Ith[j])
				CSj = 0.0;
			else
			{
				if(Ei>X[k+1])//while(X[k+1]<Ei)
					k += 1;
				if(Ei>=X[K-1])
					k = K-1;

				CSj = A[k]*Ei+B[k];
				if(CSj<0)
				{
					//printf("WARNING!!!Cross section #%d is negative at point E=%.2lfeV\n",j,Ei);
					//fprintf(log,"WARNING!!!Cross section #%d is negative at point E=%.2lfeV\n",j,Ei);
					//Serr += 1;
					CSj = 0.0;
				}
			}

			CS[i][j] = CSj;///сечения,прочитанных из файла процессов, по точкам в центре ячеек сетки

			if(add!=0)
			{
                int I;
                for(d=1;d<=add;d++)
                {
                    if(Ei<Ith[j+d])
                        CS[i][j+d] = 0.0;
                    else
                    {
                        I = i - int((Ith[j]-Ith[j+d])/dEev);//сдвиг порога по энергии
                        if(I<0)
                        {
                            fprintf(log,"\nWARNING!!!Threshold shift is larger than Ith(target) - %s Cross section #%d \n",KeyW[Kw],j+d);
                            break;
                        }
                        CS[I][j+d] = CS[i][j];

                        if(i==NEmax-1)
                        {
                            while(I<=NEmax-1)
                            {
                                CS[I][j+d] = CS[i][j];//
                                I += 1;
                            }
                        }
                    }
                }

			}

		}

		///for_reverse_process***************************************
        if(rev!=0)
        {
            int nr_t = CSR[j][3]-1;
            int nr_p = CSR[j][2]-1;
            int Ii = int(Ith[j]/dEev-0.5);

            ///Stat=g_gr/g_ex=st_w[nr_p]/st_w[nr_t];
            for(i=0;i<NEmax;i++)
            {
                Ei = (i+0.5)*dEev;

                if(i+Ii<NEmax)
                    CS[i][j+1] = CS[i+Ii][j]*(1+Ith[j]/Ei)*Stat;
                else
                    CS[i][j+1] = CS[NEmax-Ii-1][j+1];
            }

            CSR[j+1][0] = 5;///type="DEEXCITATION"
            CSR[j+1][1] = 1;///
            CSR[j+1][2] = nr_t+1;///target
            CSR[j+1][3] = nr_p+1;///product
            Ith[j+1] = Ith[j];///threshold

            ///Поправка в массив ссылок
            s = CSref[nr_t][5][0]+1;
			CSref[nr_t][5][s] = j+2;
			CSref[nr_t][5][0] += 1;
        }

		///Loging*****************************************************
		/*fprintf(log,"CS#%d at calc-energy net\n",j);
		for(i=0; i<NEmax; i++)
		{
			fprintf(log,"%.2lf\t%.2e\n",dEev*(i+0.5),CS[i][j]);
			i += 50;
		}
		fprintf(log,"\n");
		fclose(log);
		*/
		///Loging*****************************************************

		//пока не понял зачем
		//CS[j][0] = 0.0;
		//if(X[1]<dE1) {CXC[as][k][0]=Y[1];}

		j += 1+add+rev;

		fscanf(cross,"%s",CSstr);
	}
	fclose(cross);

	//Printing diagnosic message*****************************************************
	printf("EEDF Cross Sections info:\n\n");
	if(Serr>0)
	{
		printf("ATTENTION!!!Number of errors = %d\n", Serr);
		fprintf(log,"ATTENTION!!!Number of errors = %d\n", Serr);
	}
	else
	{
		printf("Number of EEDF-processes = %d\n", j);
		fprintf(log,"Number of EEDF-processes = %d\n", j);
	}

	J = j;

	//Запись в Лог***********************************************************

	fprintf(log,"Log for CS-matrix\n\n");

	fprintf(log,"E,eV\t\t");
	for(j=0;j<J;j++)
		fprintf(log,"CSnum=%d\t\t",j+1);
	fprintf(log,"\n");

	for(k=0;k<NEmax;k++)
	{
		fprintf(log,"%.2lf\t",dEev*(k+0.5));
		for(j=0;j<J;j++)
			fprintf(log,"%.2e\t",CS[k][j]);
		fprintf(log,"\n");
		//k += 10;
	}
	fprintf(log,"\n\n");

	fprintf(log,"Log for CSR[][]-matrix\n\n");

	fprintf(log,"CSnum\t");
	for(i=0;i<5;i++)
		fprintf(log,"Ind=%d\t",i);
	fprintf(log,"\n\n");

	for(j=0;j<J;j++)
	{
		fprintf(log,"%d\t",j+1);
		for(i=0;i<5;i++)
			fprintf(log,"%d\t",CSR[j][i]);
		fprintf(log,"\n");
	}
	fprintf(log,"\n\n");

	int Nt = 12;
	fprintf(log,"Log for CSref[%s][][]-matrix\n\n",Spec[Nt]);

	fprintf(log,"CStype\t\t RSum for %s\t\t CS-num\n\n",Spec[Nt]);

	int jj;
	for(Kw=0;Kw<CStype;Kw++)
	{
		jj = CSref[Nt][Kw][0];//CSref[n_t][Kw][0];
		fprintf(log,"%s\t\t%d\t",KeyW[Kw],jj);
		if(jj>0)
		{
			for(j=1;j<=jj;j++)
				fprintf(log,"%d\t",CSref[Nt][Kw][j]);
		}
		fprintf(log,"\n");
	}
	fprintf(log,"\n\n");

	//**********добавление EEDF-реакций в файл кин.схемы**********************//

	int n0,n1,n2,n3;
	FILE *react;
	react = fopen("ReactionSet.txt", "w");

	fprintf(react,"%s\n",symb);
	fprintf(react,"\nEEDF_Reactions\n");
	fprintf(react,"\n%s\n",symb);
	jj = 0;
	for(j=0;j<J;j++)
	{
		n0 = CSR[j][0];//CS-type
		n1 = CSR[j][2]-1;//target
		n2 = CSR[j][3]-1;//product-1
		n3 = CSR[j][4]-1;//product-2

		if(n0==6)//ion
		{
			fprintf(react,"R%d\t e\t + %s\t ->\t  %s\t + e\t + e\t ;\tEEDF\t//see_CS-set\n",jj+1,Spec[n1],Spec[n2]);
			jj += 1;
		}
		else if(n0==7 || n0==8)//att+rec
		{
			fprintf(react,"R%d\t e\t + %s\t ->\t %s\t ;\tEEDF\t//see_CS-set\n",jj+1,Spec[n1],Spec[n2]);
			jj += 1;
		}
		else if(n0>1)
		{
			if(n3>0)
			{
                if(n0==9)//dissatt
                {
                    fprintf(react,"R%d\t e\t + %s\t ->\t %s\t + %s ;\tEEDF\t//see_CS-set\n",jj+1,Spec[n1],Spec[n2],Spec[n3]);
                    jj += 1;
                }
                else
                {
                    fprintf(react,"R%d\t e\t + %s\t ->\t e\t + %s\t + %s ;\tEEDF\t//see_CS-set\n",jj+1,Spec[n1],Spec[n2],Spec[n3]);
                    jj += 1;
                }
			}
			else if(n2>0 && n2!=n1)
			{
				fprintf(react,"R%d\t e\t + %s\t ->\t e\t + %s ;\tEEDF\t//see_CS-set\n",jj+1,Spec[n1],Spec[n2]);
				jj += 1;
			}
		}
	}
	//fprintf(react,"\n");
	//fprintf(react,"// %s\n",symb);

	fclose(react);

	CSout = J-jj;///Elast+Rott+V-to-Vexcit_are_excluded!!!

	printf("%d of %d EEDF-processes were added to ReactionSet.txt\n\n",jj,J);
	fprintf(log,"%d of %d EEDF-processes were added to ReactionSet.txt\n\n",jj,J);

	fclose(log);

	return jj;
}
void EEDF_calc(int N,double *Ne,double *Ni,double *Te,double E,double Tgas,double Nel,double *De,double *Mue,double *DE,double *MuE,double dte,double tic,int dot,double eps)//решение уравнения Больцмана
{
	int k,n,m,s,j,J,Jmax,nte;
    double Te0,Te1,Norm,E_kT;
    double Ee,Vdr,De_Muel;

	double A,B,C,F,den;
	double Vm[NEmax],Vmi[NEmax],CSNm[NEmax],D[NEmax],Vm_av,Del,Muel,DEel,MuEel;
	double Ur,Ul,Dr,Dl,Vr,Vl;

	//Сетка по энергии:
	/*
	 Ne(k)[0]	[1]   [2]	[3]		   		 [k-1]  [k]  [k+1]						  [NEmax-1]
		|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|--x--|---------->E
	  k=0     1     2     3     4		   k-1    k    k+1                    			   NEmax
	*/

	nte = 0;
	Te1 = *Te;
	do
    //(nte=0;nte<20;nte++)
    {
        Te0 = Te1;

        //Difinition of D=A and Vm (Raizer)
        for(k=0;k<=NEmax-1;k++)
        {
            Vm[k] = 0.0;
            Vmi[k] = 0.0;
            CSNm[k] = 0.0;
            double V = 0.0;
            for(n=1;n<N;n++)//по всем компонентам смеси (с кот. сталк. эл-ны), кроме электронов
            {
                //elastic collisions
                Jmax = CSref[n][0][0];//m=0

                if(Jmax>0 && Ni[n]>0.0)
                {
                    for(j=1;j<=Jmax;j++)
                    {
                        J = CSref[n][0][j]-1;

                        //CS[k][J] = 1e-15;//Test_for_Druvestein
                        //V = 2.0e11;//Test_for_Maxwell

                        CSNm[k] += CS[k][J]*Ni[n];

                        V = CS[k][J]*sqrt(2.0/me*dE*(k+0.5))*Ni[n];//sqrt(2/me)=4.69e13

                        Vm[k] += V;//[1/s]
                        Vmi[k] += V/Mi[n];
                    }
                }
            }

            D[k] = 2*e*e*E*E/(3*me*Vm[k]);// 2*e*e/(3*me*E0*E0)=1875.46*E^2/Vm - [B/см]^2*с
            Vmi[k] = 2*me*Vmi[k];//Sum(2*me*Vi/Mi);
        }

        ///Цикл прогонки******************************************************
        ///определение прогоночных коэффициентов
        double al[NEmax],bet[NEmax];
        double Qinel,Qex;

        for(k=0;k<=NEmax-1;k++)
        {
            ///набор энергии от поля и упругие потери в соударениях
            if(k==0)
                Dl = D[k];
            else
                Dl = 0.5*(D[k]+D[k-1]);

            if(k==NEmax-1)
                Dr = D[k];
            else
                Dr = 0.5*(D[k+1]+D[k]);

            Vr = 0.5*(Vmi[k+1]+Vmi[k]);
            Vl = 0.5*(Vmi[k]+Vmi[k-1]);

            Ur = 0.5*Dr-Vr*(k+1)*dE;//point Ne[k]-right face
            Ul = 0.5*Dl-Vl*k*dE;//point Ne[k-1]-left face

            //without drift
            A = - Dl*k*dE/(dE*dE);//k-1(left)
            B = - Dr*(k+1)*dE/(dE*dE);//k+1(right)
            C = -(A + B) + 1/dte;//k(center)

            //with drift part
            if((Ur>=0)&&(Ul>=0))
            {
                A += -Ul/dE;
                C += Ur/dE;
            }
            if((Ur>=0)&&(Ul<0))
                C += (Ur-Ul)/dE;
            if((Ur<0)&&(Ul<0))
            {
                C += -Ul/dE;
                B += Ur/dE;
            }
            if((Ur<0)&&(Ul>=0))
            {
                A += -Ul/dE;
                B += Ur/dE;
            }

            ///inelastic collisions*************************************************************************************
            Qinel = 0.0;
            s = 0;
            for(n=1;n<N;n++)//суммирование по всем компонентам смеси (с кот. сталк. эл-ны)
            {
                if(Ni[n]>0.0)
                {
                    for(m=2;m<CStype;m++)//суммирование по всем типам процессов
                    {
                        Jmax = CSref[n][m][0];//число пар (Ntarget,CStype)

                        if(Jmax>0)//только по ненулевым типам процессов
                        {
                            if((m>=2)&&(m<=4))///excitation(rot,vib,elec,diss)
                            {
                                for(j=1;j<=Jmax;j++)
                                {
                                    J = CSref[n][m][j]-1;//номер сечения для пар (Ntarget,CStype)

                                    s = k+int(Ith[J]/dEev);
                                    if(s>=NEmax-1)
                                        s = NEmax-1;

                                    Qex = -CS[k][J]*sqrt(2/me*dE*(k+0.5))*Ne[k]+CS[s][J]*sqrt(2/me*dE*(s+0.5))*Ne[s];

                                    Qinel += Qex*Ni[n];
                                }
                            }
                            if(m==5)///deexcitation
                            {
                                for(j=1;j<=Jmax;j++)
                                {
                                    J = CSref[n][m][j]-1;

                                    s = k-int(Ith[J]/dEev);

                                    Qex = -CS[k][J]*sqrt(2/me*dE*(k+0.5))*Ne[k];
                                    if(s>0)
                                        Qex += CS[s][J]*sqrt(2/me*dE*(s+0.5))*Ne[s];
                                    Qinel += Qex*Ni[n];
                                }
                            }
                            if(m==6)///ionization_2-models
                            {
                                for(j=1;j<=Jmax;j++)
                                {
                                    J = CSref[n][m][j]-1;

                                    /*//равнораспределение между 2-мя эл-нами********************************************
                                    s = 2*k+int(Ith[J]/dEev);
                                    if(s>=NEmax-1)
                                        s = NEmax-1;

                                    Qex = -CS[k][J]*sqrt(2/me*dE*(k+0.5))*Ne[k]+2*CS[s][J]*sqrt(2/me*dE*(s+0.5))*Ne[s];
                                    */
                                    //*********************************************************************************

                                    //ионизирующий эл-н уносит остаток, новый уносит E0-Райзер+Хагеллар****************
                                    E_kT = Tgas/eV_K;//E_kT=kT(Волошин)[эВ]//
                                    s = k+int((Ith[J]+E_kT)/dEev);
                                    if(s>=NEmax-1)
                                        s = NEmax-1;

                                    Qex = 0.0;
                                    if((k>=int(E_kT/dEev-0.5))&&((k<=int(E_kT/dEev+0.5))))
                                    {
                                        int kI = int((Ith[J]+E_kT)/dEev);
                                        int ki = 0;
                                        for(ki=kI;ki<NEmax;ki++)
                                            Qex += CS[ki][J]*sqrt(2/me*dE*(ki+0.5))*Ne[ki];
                                    }

                                    Qex = -CS[k][J]*sqrt(2/me*dE*(k+0.5))*Ne[k]+CS[s][J]*sqrt(2/me*dE*(s+0.5))*Ne[s]+Qex;
                                    //**********************************************************************************

                                    Qinel += Qex*Ni[n];
                                }
                            }
                            if(m>=7)///attachment+recomb+diss-attach
                            {
                                for(j=1;j<=Jmax;j++)
                                {
                                    J = CSref[n][m][j]-1;

                                    Qex = -CS[k][J]*sqrt(2/me*dE*(k+0.5))*Ne[k];
                                    Qinel += Qex*Ni[n];
                                }
                            }


                            //rotation??
                        }
                    }
                }
            }
            ///*********************************************************************************************************

            ///RH-part
            F = Ne[k]/dte + Qinel;

            if(k==0)
            {
                al[k+1] = -B/C;
                bet[k+1] = F/C;
            }
            else if(k==NEmax-1)
            {}
            else
            {
                den = 1.0/(A*al[k]+C);
                al[k+1] = -B*den;
                bet[k+1] = (F-A*bet[k])*den;
            }
        }

        ///граничное условие**************************************************
        Ne[NEmax-1] = 0.0;
        //Ne[NEmax-1] = (F-A*bet[NEmax-1])/(A*al[NEmax-1]+C);

        ///цикл обратной прогонки*********************************************
        for(k=NEmax-2;k>=0;k--)
        {
            Ne[k] = al[k+1]*Ne[k+1]+bet[k+1];
            if(Ne[k]<1e-30)
                Ne[k] = 0.0;
        }

        ///Normalization***********************************************
        Norm = 0.0;
        for(k=0; k<NEmax; k++)
            Norm += Ne[k]*dEev;

        for(k=0; k<NEmax; k++)
            Ne[k] = Ne[k]*Nel/Norm;

        ///Ee_Te-calculation*******************************************
        Ee = 0.0;
        for(k=0;k<NEmax;k++)
            Ee += Ne[k]*(k+0.5)*dEev*dEev;
        Ee = Ee/Nel;//[eV]
        Te1 = 2.0/3.0*Ee;//[eV]

        nte++;

    }while(fabs(1.0-(Te1/Te0))>eps);// || (nte<100));//(fabs(Te1-Te0)>0.01);////(fabs(Te1-Te0)>0.001);

    ///Apply_to_used_matrices******************************************
    *Te = Te1;

    ///De_&_Mue-calculation********************************************
    Del = 0.0;
    for(k=0;k<NEmax;k++)
        Del += sqrt((k+0.5)*dE)*Ne[k]/CSNm[k];
    Del *= sqrt(2.0/me)/(3.0*Nel)*dEev;
    *De = Del;///[см2/с]

    Muel = 0.0;
    for(k=0;k<NEmax-1;k++)
        Muel += (k+1)*(Ne[k+1]/sqrt(k+1.5)-Ne[k]/sqrt(k+0.5))/CSNm[k];
    Muel *= -sqrt(2.0*dE/me)*e/(3.0*Nel)*dEev/dE;
    *Mue = Muel;///[cm2/c*abs]

    De_Muel = Del*Eabs/Muel;//[В]
    Vdr = Muel*E;///[см/с]
    //*Jel = e*Nel*Vdr;//[Abs/cm2*s]

    //Qel = me/Mi[11]*Ee*Vm_av*Nel*1.602e-12;//[erg/cm3*s]

    //QE = E*Jel*3.14*pi*Rad*Rad; //[abs*СГС/s]

    ///DE_&_MuE-calculation********************************************
    DEel = 0.0;
    for(k=0;k<NEmax;k++)
        DEel += pow((k+0.5)*dE,1.5)*Ne[k]/CSNm[k];
    DEel *= sqrt(2.0/me)/(3.0*Nel*Ee)*dEev*dEev/dE/Ee;
    *DE = DEel;///[см2/с]

    MuEel = 0.0;
    for(k=0;k<NEmax-1;k++)
        MuEel += (k+1)*(k+1)*(Ne[k+1]/sqrt(k+1.5)-Ne[k]/sqrt(k+0.5))/CSNm[k];
    MuEel *= -sqrt(2.0*dE/me)*e/(3.0*Nel*Ee)*dEev*dEev/dE/Ee;
    *MuE = MuEel;///[cm2/c*abs]


    ///Writing_data***************************************************
	//if(dot==Ndots)
		//EEDF_print(Ne,*Te,Nel,Norm,tic);

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

}
void EEDF_const_calc(char *Txt,int Nedf,double *Kel,double Nel,double E,double Ngas,double *Te,double *dTe,double *De,double *Mue,double tic)//Kel(EEDF)-calculation
{
	int n,m,s,Smax,sL,sR;
    double Te0,Ki,Y[Nedf+5],dEN,A,B,EN,ddTe;

    Te0 = *Te;
    EN = fabs(E)/Ngas*1e17*Eabs;
    Smax = int(EN_Tbl/dEN_Tbl);

    ///From_EEDF_Table***********************************
    if(!strcmp(Txt,"E/N-case"))///E/N-case***************
    {
        if(EN>EN_Tbl)
        {
            for(m=0;m<Nedf+7;m++)
                Y[m] = EedfTbl[Smax][m+2];
        }
        else if(EN<EedfTbl[0][m+2])
        {
            for(m=0;m<Nedf+7;m++)
                Y[m] = EedfTbl[0][m+2];
        }
        else
        {
            s = int(EN/dEN_Tbl);
            if(EN>EedfTbl[s][0] && EN<EedfTbl[s+1][0])
                s++;

            dEN = dEN_Tbl;//EedfTbl[s][0]-EedfTbl[s-1][0];
            for(m=0;m<Nedf+7;m++)
            {
                A = (EedfTbl[s+1][m+2]-EedfTbl[s][m+2])/dEN;
                B = EedfTbl[s][m+2]-A*EedfTbl[s][0];
                Y[m] = A*EN+B;
            }

            /*
            if(s==0)
            {
                for(m=0;m<Nedf+7;m++)
                    Y[m] = EedfTbl[0][m+2];
            }
            else if(s>0 && EN<=EN_Tbl)
            {
                dEN = dEN_Tbl;//EedfTbl[s][0]-EedfTbl[s-1][0];
                for(m=0;m<Nedf+7;m++)
                {
                    A = (EedfTbl[s][m+2]-EedfTbl[s-1][m+2])/dEN;
                    B = EedfTbl[s-1][m+2]-A*EedfTbl[s-1][0];
                    Y[m] = A*EN+B;
                }
            }
            */
        }

        //*Te = Y[0];
        *dTe = fabs(*Te-Te0);
    }
    else if(!strcmp(Txt,"Te-case"))///Te-case**************!!!!!!!!CHECKKKKKKKKKKKKKKKKKKKKKKKKK!!!!!!!!!!!!!!!!!!
    {
        if(Te0>EedfTbl[Smax][2])
        {
            for(m=0;m<Nedf+7;m++)
                Y[m] = EedfTbl[Smax][m+2];
        }
        else if(Te0<EedfTbl[0][2])
        {
            for(m=0;m<Nedf+7;m++)
                Y[m] = EedfTbl[0][m+2];
        }
        else
        {
            sL = 0,sR = Smax;
            do
            {
                s = int((sL+sR)/2);

                if(Te0>EedfTbl[s][2])
                    sL = s;
                else
                    sR = s;
            }while(!((Te0>EedfTbl[s][2]) && (Te0<EedfTbl[s+1][2])));

            ddTe = EedfTbl[s+1][2]-EedfTbl[s][2];
            for(m=0;m<Nedf+7;m++)
            {
                A = (EedfTbl[s+1][m+2]-EedfTbl[s][m+2])/ddTe;
                B = EedfTbl[s][m+2]-A*EedfTbl[s][0];
                Y[m] = A*Te0+B;
            }
        }
    }

    ///Присвоение_значений используемым переменным
    *Mue = Y[1]/Ngas*Eabs;//Mue
    *De = Y[2]/Ngas;//1.5;//Proshina_check!
    *(Mue+(LEN+2)) = Y[3]/Ngas*Eabs;///MueE
    *(De+(LEN+2)) = Y[4]/Ngas;///DeE
    for(m=0;m<Nedf;m++)
        Kel[m] = Y[m+5];
    //Kel[16] /= 1.5;////Proshina_check!

    /*
    ///Logging*******************************************
    FILE *log;
    if(tic==0)
        log = fopen("Kel_data.txt", "w");
    else
        log = fopen("Kel_data.txt", "a+");

    fprintf(log,"t=%.2e[s]\n",tic);
	for(m=0; m<Nedf; m++)
	    fprintf(log,"%d\t%.2e\n",m+1,Kel[m]);
	fprintf(log,"\n");
	fclose(log);
	///Logging*******************************************
	*/
}
void EEDF_print(double *Ne,double Te,double Nel,double Norm,double tic)//запись EEDF в файл
{
	int k;

    /*
	printf("\nt = %.2e[s]\n",tic);
	printf("Te = %.2lf[eV]\t Ne = %.2e[cm-3]\n",Te,Nel);
	*/

	FILE *log;
	log = fopen("EEDF_data.txt", "w");

    fprintf(log,"t=%.2e[s]\t Te=%.2lf[eV]\n",tic,Te);
    fprintf(log,"Ne=%.2e[cm-3]\n",Nel);
	for(k=0; k<NEmax; k+=10)
        fprintf(log,"%.2lf\t%.2e\n",(k+0.5)*dEev,Ne[k]/Nel);
	fprintf(log,"\n\n");

	/*
	fprintf(log,"t=%.2e [s]\nNorm=%.2e [cm-3]\nTe=%.2lf [eV]\nE,[eV]\t Ne(E),[cm-3/eV]\n",dt*(nt+1),Norm,Te);//Nel
	for(k=0; k<NEmax; k+10)
		fprintf(log,"%.2lf\t%.2e\n",dEev*(k+0.5),Ne[k]/Nel);
	fprintf(log,"\n\n");
	*/
	fclose(log);
}
void EEDF_table_calc(int Nedf,int N,int NLen,double *Nn,double Ng,double Tg,double eps,char *EedfTblFl)
{
    double E,E_N,Te,Nel,Del,Muel,DEel,MuEel,E_Nmin,Vdr,Ki;
    double Kel[Nedf],Ni[N],Ne[NEmax];
    int st,m,n,k,I,dot = 0;

    E_Nmin = 0.5;///[Td]

    for(n=0;n<N;n++)
        Ni[n] = *(Nn+n*NLen);
    Nel = Ni[0];

    ///Initial_EEDF****************************************************
    Te = 0.5;
    for(k=0;k<NEmax;k++)
    {
        Ne[k] = 2*Nel*pow((k+0.5)*dEev/(Te*Te*Te*pi),0.5)*exp(-(k+0.5)*dEev/Te);//Maxwell
        if(Ne[k]<1e-30)
            Ne[k] = 0.0;
    }
	///****************************************************************

    EN_Tbl = 1000.0;
    dEN_Tbl = 5.0;
    E_N = E_Nmin;
    I = 0;
    st = 0;
    do
    {
        E = E_N*Ng*1e-17/Eabs;

        EEDF_calc(N,Ne,Ni,&Te,E,Tg,Nel,&Del,&Muel,&DEel,&MuEel,1.0e-11,0.0,dot,eps);

        for(m=0;m<Nedf;m++)
        {
            Ki = 0.0;
            for(k=0;k<NEmax;k++)//int(Ith[m])
                Ki += CS[k][m+CSout]*sqrt(2.0/me*dE*(k+0.5))*Ne[k]*dEev;
            Kel[m] = Ki/Nel;
        }

        //EEDF_print(Ne,Te,Nel,E_N);

        Vdr = Muel*E_N*Ng*1e-17/Eabs;

        ///EedfTable_Data_Saving*************************************************
        EedfTbl[I][0] = E_N;
        EedfTbl[I][1] = Vdr;
        EedfTbl[I][2] = Te;
        EedfTbl[I][3] = Muel*Ng/Eabs;//[1/(cm*s*V)]
        EedfTbl[I][4] = Del*Ng;//[1/(cm*s)]
        EedfTbl[I][5] = MuEel*Ng/Eabs;//[1/(cm*s*V)]
        EedfTbl[I][6] = DEel*Ng;//[1/(cm*s)]
        for(int m=0;m<Nedf;m++)
            EedfTbl[I][m+7] = Kel[m];

        if(I==0)
            E_N = 0.0;
        E_N += dEN_Tbl;

        if(E_N>EN_Tbl && st!=0)
        {
            E_N = EN_Tbl;
            st++;
        }

        I++;

    }while(E_N<=EN_Tbl);

    EEDF_table_print(Nedf,I,EedfTblFl);
}
void EEDF_table_print(int Nedf,int I,char *EedfTblFl)
{
    FILE *tbl;
    int i,m;

    tbl = fopen(EedfTblFl, "w");
    fprintf(tbl,"E/N,Td\tVdr,cm/s\tTe,eV\tMue*N,1/(cm*s*V)\tDe*N,1/(cm*s)\t\tMuE*N,1/(cm*s*V)\tDE*N,1/(cm*s)\t");
    for(m=0;m<Nedf;m++)
        fprintf(tbl,"%s\t",RName[m]);
    fprintf(tbl,"\n");

    for(i=0;i<I;i++)
    {
        fprintf(tbl,"%.2lf\t",EedfTbl[i][0]);
        for(m=1;m<Nedf+7;m++)
            fprintf(tbl,"%.3e\t",EedfTbl[i][m]);
        fprintf(tbl,"\n");
    }

    printf("EEDF_Table_has_been_created_File:O2_EEDF_Ke_Table.txt\n\n");

    fclose(tbl);
}
int EEDF_table_read(int N,int Nedf,char *EedfTableFile)
{
    FILE *tbl;
    tbl = fopen(EedfTableFile, "r");

    FILE *log;
	log = fopen("Log_EedfTable.txt", "w");

    char Rstr[50];
    int i,n,err;

    if(!tbl)
    {
        err++;

        printf("\nNO EEDF-data file with %s name!!\n",EedfTableFile);
        printf("EEDF-data will be recalculated!\n");

        fprintf(log,"NO EEDF-data file with %s name!!\n",EedfTableFile);
        fprintf(log,"EEDF-data will be recalculated!\n");
    }
    else
    {
        fprintf(log,"The EEDF-file %s has been scanned\n",EedfTableFile);
        ///Scanning until the name of processes********************************************
        do
        {
            fscanf(tbl,"%s",Rstr);
        }while(strcmp("DE*N,1/(cm*s)",Rstr));

        ///Scanning names of processes and check them on matching with CS-data*************
        err = 0;
        for(n=0;n<Nedf;n++)
        {
            fscanf(tbl,"%s",Rstr);
            if(strcmp(Rstr,RName[n]))
            {
                fprintf(log,"Reaction R#%d - %s is not matching the e-Process %s from CS-file\n",n+1,Rstr,RName[n]);
                err++;
            }
        }

        ///Saving_the_data_from_EEDF-Table*************************************************
        if(err==0)
        {
            i = 0;
            do
            {
                for(n=0;n<Nedf+7;n++)
                    fscanf(tbl,"%lf",&EedfTbl[i][n]);

                EN_Tbl = fmax(EN_Tbl,EedfTbl[i][0]);

                i++;
            }while(i<300);
            dEN_Tbl = EedfTbl[2][0] - EedfTbl[1][0];

            printf("\nThe EEDF-data has been Successfully imported from %s!!!\n\n",EedfTableFile);
            fprintf(log,"The EEDF-data has been Successfully imported!!!\n");
        }
        else
        {
            printf("\nWarning!!!Problems while EEDF-data import have been detected!See Log_EedfTable.txt\n");
            printf("EEDF-data will be recalculated!\n\n");
            fprintf(log,"EEDF-data will be recalculated!\n");
        }
    }

	fclose(tbl);
	fclose(log);

	return err;
}
