#define CollisionalFreq 0//

double PlasmaFreq(int z)
{
return 10e6*2.0*pii;//1.0e7*2.0*pii;//sqrt(5.0)*17.83e9/2;sqrt((z-MZ/2+thick)/(2.0*thick));//wp.
}

double chi(int m,int z)
{
	double wp,vc;
	wp=PlasmaFreq(z);
	vc=CollisionalFreq;
//	return wp*wp*dt/vc-(wp*wp/(vc*vc))*(1.0-exp(-vc*dt))*exp(-m*vc*dt);
	return (epsn*wp*wp/vc)*(1.0-exp(-1.0*vc*dt))*exp(-1.0*m*vc*dt);

}

double xi(int m,int z)
{
	double wp,vc;
	wp=PlasmaFreq(z);
	vc=CollisionalFreq;
	return (epsn*wp*wp/(vc*vc*dt))*(1.0-(1.0+vc*dt)*exp(-1.0*vc*dt))*exp(-m*vc*dt);

 /*   double wp,vc,xi;
	wp=PlasmaFreq(z);
	vc=CollisionalFreq;
	xi=wp*wp*dt/(2.0*vc);
	return xi-wp*wp/(vc*vc*vc*dt)*(1.0-(1.0+vc*dt)*exp(-vc*dt))*exp(-m*vc*dt);
*/
	//return 0.0;
}

inline double dChi0(int z){	return chi(0,z)-chi(1,z);}
inline double dXi0(int z){	return xi(0,z)-xi(1,z);}
int tri_matrix(int n, double *a, double *b, double *c, double *f)
{
	int kk;
	double w;
	f[0]=f[0]/b[0];
	w=b[0];
	for (kk=1; kk<=n-1;kk++)
	{
		b[kk-1]=c[kk-1]/w;
		w=b[kk]-a[kk-1]*b[kk-1];
		f[kk]=(f[kk]-a[kk-1]*f[kk-1])/w;
	}
	for(kk=n-1; kk>=1; kk--) f[kk-1]=f[kk-1]-b[kk-1]*f[kk];
	return(1);
}