#include "peg_pic.h"

float *gam;
float source();       /* calculates U_bias, or I_total, 
                    depending on the circuit type (v or i)*/
void  rpot();	   // Reads the specified potential (if Field=2)

/********************************************************/
/* Solve a tridiagnal matrix to get Phi[i] at time t.   */
/* bcptr=d_0 (or d_1 for short circuit), see ref. below.*/
/* circuit = sigma (open cirquit), = I (short circuit), */ 
/* or calculates sigma and total and circuit currents   */
/*                             for the general case.    */

void fields() {
	static int init_flag=1;
	int isp, i;
	float bet, s; 

	if(init_flag) {
		if ( field == 2 ) {
			/* read potential from file */
			rpot();
		} else {
			gam  = (float *)malloc(ng*sizeof(float));
		}

		init_flag = 0;
	}

	jp = backj;
	
	for (isp=0; isp<nsp; isp++) {
		jp += jwall[isp];
	}
	
	
	if (field < 2 || field == 4 ) {
		/* Calculates the self-consistent potential
		 * in case of self-consistent field or first step when using phi from previous step */
		if (nstrt) {
			phi[0] = source();
		}
		
		bet = b0;
		phi[nstrt] = (*bcptr)()/bet;  

		if (cappa) {
			for (i=nstrt+1; i<nc; i++) {
				phi[i] = se*rho[i] +phi[i-1];
			}

			phi[nc] = 0.;
			
			for(i=nc-2; i>=nstrt; i--) {
				phi[i] += phi[i+1];
			}
		} else {
			for (i=nstrt+1; i<nc; i++) {
				gam[i] = 1./bet;
				bet = -2 - gam[i];
				phi[i] = -(se*rho[i] +phi[i-1])/bet;
			}

			phi[nc] = 0.;
			
			for(i=nc-2; i>=nstrt; i--) {
				phi[i] -= gam[i+1]*phi[i+1];
			}
		}

		(*circuitptr)();
		// end of if (field<2) routine
	} 
	
	/**** Calculate E(x) from phi(x) ****/

	s = 0.5/dx;
	
	for (i=1; i<nc; i++) {
		e[i] = (phi[i-1] - phi[i+1])*s;
	}
	
	e[0]  = (phi[0] - phi[1])/dx - rho[0]*0.5*dx/epsilon;
	e[nc] = (phi[nc-1] - phi[nc])/dx + rho[nc]*0.5*dx/epsilon;

}   /* end of FIELDS */

/***************************************************************/
/***              For field solver                           ***/
/** see J.P. Verboncouer, J. Comp. Phys., 104 (2), 321 (1993) **/

float source() {
	if (dcramped) {
		if(t < risetime) {
			if(ramp) {
				return (ramp*t);
			} else {
				return (.5*dcbias*(1- cos(t*w0)));
			}
		} else {
			return (dcbias);
		}
	} else {
		return (dcbias + ramp*t + acbias*sin(t*w0 + theta0));
	}
}


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

/* When C = 0: OPEN CIRCUIT */

float bc1() {
	return (-(sigma + dt*jp + 0.5*rho[0]*dx)*dx/epsilon);
}

void circuit1() {
	oldsigma = sigma;
	sigma += jp*dt;
}

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

/* When CURRENT SOURCE is applied */

float bc2() {
	return (-(sigma + dt*(source()/area + jp) + 0.5*rho[0]*dx)*dx/epsilon);
}

void circuit2() {
	exti  = source();
	oldsigma = sigma;
	sigma+= dt*(jp + exti/area);
}

/***************************************************************/
/* When R=L=0, C -> infinity: SHORT CIRCUIT */

float bc3() {
	return (-rho[1]*dx*dx/epsilon - source());
}

void circuit3() {
	oldsigma = sigma;
	sigma= (phi[0] - phi[1])*epsilon/dx - rho[0]*dx*0.5;
	exti = area*((sigma - oldsigma)/dt - jp);
}

/***************************************************************/
/* The General case */

float bc4() {
	static float a0, a1, a2, a3, a4;
	float k;

	if(!a0) {
		a0 = 2.25*extl/dt/dt + 1.5*extr/dt + 1/extc;
		a1 = -6*extl/dt/dt - 2*extr/dt;
		a2 = 5.5*extl/dt/dt + .5*extr/dt;
		a3 = -2*extl/dt/dt;
		a4 = .25*extl/dt/dt;
	}
	
	k = (a1*extq + a2*extq_1 + a3*extq_2 + a4*extq_3)/a0;
	
	return (-(0.5*rho[0]*dx + sigma + dt*jp
			+ (source()/a0 - k - extq)/area)*dx/epsilon);
}

void circuit4() {
	static float a0, a1, a2, a3, a4;
	float k;

	if(!a0) {
		a0 = 2.25*extl/dt/dt + 1.5*extr/dt + 1/extc;
		a1 = -6*extl/dt/dt - 2*extr/dt;
		a2 = 5.5*extl/dt/dt + .5*extr/dt;
		a3 = -2*extl/dt/dt;
		a4 = .25*extl/dt/dt;
	}
	
	k = (a1*extq + a2*extq_1 + a3*extq_2 + a4*extq_3)/a0;
	extq_3 = extq_2;
	extq_2 = extq_1;
	extq_1 = extq;
	extq = (source() - phi[0])/a0 - k;

	exti = (extq - extq_1)/dt;
	oldsigma= sigma;
	sigma  += (jp + exti/area)*dt;
}

/***************************************************************/
/* Reading the file with the potential */

void rpot() {
	int            i;
	float          ttmp;
	FILE           *mytxt;

	if ((mytxt = fopen(dargv5, "r")) == NULL) {
		puts("The file potential is missing!");
		exit(1);
	}

	for(i=0; i<=nc; i++) {
		fscanf(mytxt, "%e\n", &ttmp);
		phi[i]= ttmp;
	}

	fclose(mytxt);

} // End of void  rpot()

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

