#include <stdlib.h>
#include <malloc.h>
#include "pdp1.h"

void freqanalysis(void);
extern void realft(float data[], int n, int isign);


// Function double frandbase2(void) {{{
double frandbase2(void)
{
	static int counter2=0;
	return (revers_base(counter2++, 2));
} // }}}


// Function double revers_base(int num, int n) {{{
double revers_base(int num, int n)
{
	double power, rev;
	int inum, iquot, irem;

	rev = 0.;
	inum = num;
	power = 1.;

	do {
		iquot = inum/n;
		irem = inum - n*iquot;
		power /= n;
		rev += irem*power;
		inum = iquot;
	} while (inum > 0);

	return (rev);
} // }}}

//**********************************************************************
// irand(iseed) returns values 1 through 2**31-2.
// From "Random number generators: good ones are hard to find", S. Park
// and K. Miller, Communications of ACM, October 1988, pp 1192-1201.
// This is from page 1195, and is to work on any system for which
// maxint is 2**31-1 or larger. Due earlier to Schrage, as cited by P&M.
//
// Note: OK values for iseed are 1 through 2**31-2. Give it 0 or 2*31-1
// and it will return the same values thereafter!
//
// C version 6/91, Bruce Langdon.
//
// Algorithm replaces seed by mod(a*seed,m). First represent
// seed = q*hi + lo.  Then
// a*seed = a*q*hi + lo = (m - r)*hi + a*lo = (a*lo - r*hi) + m*hi,
// and new seed = a*lo - r*hi unless negative; if so, then add m.

// Function float frand(void) {{{
float frand(void)
{
	long a = 16807, m = 2147483647, q = 127773, r = 2836;
	long hi, lo;
	float fnumb;

	hi = seed/q;
	lo = seed - q*hi;
	seed = a*lo - r*hi;
	// "seed" will always be a legal integer of 32 bits (including sign).
	if (seed <= 0) seed = seed + m;
	fnumb = seed/2147483646.0;

	return (fnumb);
} // }}}

// Returns bit reversed num in base 2

// Function double revers(unsigned int num) {{{
double revers(unsigned int num)
{
	double f=0.5, sum=0.0;

	while(num) {
		if (num & 1) sum += f;     // is 1st bit set?
		f *= 0.5;
		num >>= 1;		       // fast divide by 2
	}
	return (sum);
} // }}}


// Function double fdist(double vv, int k, int isp) {{{
double fdist(double vv, int k, int isp)
{
	double u;

	u = (vv - v0[isp][k])/vt[isp][k];
	return (exp(-u*u*0.5));
} //  }}}

//**********************************************************************
// time history accumulator; calculates and stores all history values,
// and performs combing on history values when low on memory

// Function void history(void) {{{
void history(void)
{
	register int i, j, isp, k;
	static float jtemp[NSMAX];
	static int count=1, nrftcount=0, nrfcount=0;

	//*****   MCC rates  *****
	if (it[0]) {
		for (i=0; i<ndiag; i++)
			for (j=0; j<ng; j++)
				rate_show[i][j] =
					(rate_show[i][j]*(it[0]-1)
					 + mccrate[i][j])/((float)it[0]);
	}

	//  Fixing the diagnostic arrays {{{
	if (nsp && it[0]) {
		for (isp=0; isp<nsp; isp++) {
			jtemp[isp] += dt*jwall[isp];

			for (i=0; i<nbin_mid[isp]; i++)
				fe_mid_show[isp][i] =
					fe_mid[isp][i] / (it[isp]*sqrt(e_mid_array[isp][i]+DBL_MIN));

			for (i=0; i< ng; i++) {
				if (sp_n_k[isp][i] <1e-10) {
					sp_u_show[isp][i] = 0.0;
					sp_ke_show[isp][i]= 0.0;
				}
				else {
					sp_u_show[isp][i] = u_scale[isp]*sp_j[isp][i]/sp_n_k[isp][i];
					sp_ke_show[isp][i]= Escale[isp]*sp_ke[isp][i]/sp_n_k[isp][i];
				}
				jdote_show[isp][i]= jdote[isp][i]*jdote_scale[isp]/it[isp];
				sp_j_show[isp][i] = sp_j[isp][i]*j_scale[isp];
			}
		}
	} // }}}

	if (nfft) {
		if (thist_hi >= nfft) {
			freqanalysis();
			thist_hi=0;
		}
		cur_hist[thist_hi]= exti;
		pow_hist[thist_hi]= exti*phi[0];
		phi_hist[0][thist_hi] = phi[0];
		phi_hist[1][thist_hi] = phi[nc/2];
		Local_t_array[thist_hi] = t;
		thist_hi++;

		// Calculating the time ave. moments {{{

		if (nrftcount == nfft) {
			nrfcount++;

			for (isp=0; isp<nsp; isp++)
				for (i=0; i<ng; i++) {
					sp_n_ave_show[isp][i]	=	sp_n_ave[isp][i]/nfft;
					sp_n_ave[isp][i]		=	sp_n[isp][i];
					sp_u_ave_show[isp][i]	=	sp_u_ave[isp][i]/nfft;
					sp_u_ave[isp][i]		=	sp_u_show[isp][i];
					sp_j_ave_show[isp][i]	=	sp_j_ave[isp][i]/nfft;
					sp_j_ave[isp][i]		=	sp_j_show[isp][i];

					sp_ke_ave_show[isp][i]	=	sp_ke_ave[isp][i]/nfft;
					sp_ke_ave[isp][i]		=	sp_ke_show[isp][i];
					sp_ke_ave_show[isp][i]	-=	0.5 * m[isp] *
						sp_u_ave_show[isp][i] *
						sp_u_ave_show[isp][i] / 1.602e-19;
					sp_ke_ave_show[isp][i] *= 2.0/3.0;
				}

			for (i=0; i<ng; i++) {
				e_ave_show[i]= e_ave[i]/nfft;      e_ave[i]  = e[i];
				phi_ave_show[i]= phi_ave[i]/nfft;  phi_ave[i]= phi[i];
			}
			nrftcount =0;
		}
		else {
			for (isp=0; isp<nsp; isp++)
				for (i=0; i<ng; i++) {
					sp_n_ave[isp][i] += sp_n[isp][i];
					sp_u_ave[isp][i] += sp_u_show[isp][i];
					sp_j_ave[isp][i] += sp_j_show[isp][i];
					sp_ke_ave[isp][i]+= sp_ke_show[isp][i];
				}
			for (i=0; i<ng; i++) {
				e_ave[i]   += e[i];
				phi_ave[i] += phi[i];
			}
			nrftcount++;
		}
	} // }}}

	if (--count) return;		// only accum every interval steps
	if (hist_hi >= HISTMAX) {	// comb time histories
		for (isp=0; isp<nsp; isp++) {
			for (i=1, k=4; i<HISTMAX/4; i++, k+=4) {
				np_hist[isp][i] = np_hist[isp][k];
				kes_hist[isp][i] = kes_hist[isp][k];
				jwall_hist[isp][i]= jwall_hist[isp][k];
			}
		}
		for (i=1, k=4; i<HISTMAX/4; i++, k+=4) {
			com_phi_hist[0][i] = com_phi_hist[0][k];
			com_phi_hist[1][i] = com_phi_hist[1][k];
			com_pow_hist[i] = com_pow_hist[k];
			com_cur_hist[i] = com_cur_hist[k];
			wall_sigma_hist[i] = wall_sigma_hist[k];
			ese_hist[i] = ese_hist[k];
			t_array[i] = t_array[k];
		}
		hist_hi = i;
		interval *= 4;
	}
	for (isp=0; isp<nsp; isp++) {		// accumulate histories
		jwall_hist[isp][hist_hi]= jtemp[isp];
		jtemp[isp]= 0.0;
		if (hist_hi)  jwall_hist[isp][hist_hi] +=jwall_hist[isp][hist_hi-1];
		np_hist[isp][hist_hi] = np[isp];
		kes_hist[isp][hist_hi] = 0;
		for (i=0; i<np[isp]; i++)
			kes_hist[isp][hist_hi] += vx[isp][i]*vx[isp][i] +vy[isp][i]*vy[isp][i]
				+vz[isp][i]*vz[isp][i];

		if (np[isp]) kes_hist[isp][hist_hi] *= Escale[isp]/np[isp];
	}
	t_array[hist_hi] = t;
	wall_sigma_hist[hist_hi] = sigma;
	com_cur_hist[hist_hi] = exti;
	com_phi_hist[0][hist_hi] = phi[0];
	com_phi_hist[1][hist_hi] = phi[nc/2];
	com_pow_hist[hist_hi] = exti*phi[0];

	ese_hist[hist_hi] = 0.5*(rho[0]*phi[0] + rho[nc]*phi[nc]);
	for (i=1; i<nc; i++) ese_hist[hist_hi] += rho[i]*phi[i];
	ese_hist[hist_hi] *= 0.5*area*dx;
	ese_hist[hist_hi] -= epsilon*e[nc]*area*phi[nc];

	ese_hist[hist_hi] = fabs(ese_hist[hist_hi]+DBL_MIN);
	for (isp=0; isp<nsp; isp++)
		kes_hist[isp][hist_hi] = kes_hist[isp][hist_hi]+DBL_MIN;

	hist_hi++;
	count = interval;
} // }}}

// Function void freqanalysis(void) {{{
void freqanalysis(void)
{
	register int i, j;
	static int init_freq_flag=1;
	static float *temp1, *temp2;

	if (init_freq_flag) {
		if (!(temp1 = (float *) malloc((nfft/2)*sizeof(float))) ||
								!(temp2 =
									(float *)malloc((nfft/2)*sizeof(float))))
			puts("Null ptr in freqanalysis()");
		init_freq_flag=0;
	}

	//****************************************************
	for (i=0; i< nfft; i++) {
		cur_fft[i]= cur_hist[i];
		pow_fft[i]= pow_hist[i];
		phi_fft[i]= phi_hist[0][i];
		mphi_fft[i]= phi_hist[1][i];
	}

	realft(phi_fft-1, freq_hi, 1);
	realft(mphi_fft-1,freq_hi, 1);
	realft(cur_fft-1, freq_hi, 1);
	realft(pow_fft-1, freq_hi, 1);

	//*** Computing mag and phase of the current signal ** {{{

	temp1[0]= fabs(cur_fft[0])/freq_hi/2;
	temp2[0]= (cur_fft[0] > 0.0) ? 0.0 : 180.0;
	for (i=1, j=2; i< freq_hi; i++, j+=2) {
		temp1[i]= sqrt(cur_fft[j]*cur_fft[j]+ cur_fft[j+1]*cur_fft[j+1])/freq_hi;
		if (fabs(cur_fft[j+1]) < 1e-30 && fabs(cur_fft[j]) < 1e-30)
			temp2[i]= 0.0;
		else
			temp2[i]= (180./M_PI)*atan2(cur_fft[j], cur_fft[j+1]);
	}
	for (i=0; i< freq_hi; i++) {
		cur_fft[i]= temp1[i];
		cur_fft[freq_hi+i]= temp2[i];
	} // }}}

	//*** Computing mag and phase of the power signal **** {{{

	temp1[0]= fabs(pow_fft[0])/freq_hi/2;
	temp2[0]= (pow_fft[0] > 0.0) ? 0.0 : 180.0;
	for (i=1, j=2; i< freq_hi; i++, j+=2) {
		temp1[i]= sqrt(pow_fft[j]*pow_fft[j]+ pow_fft[j+1]*pow_fft[j+1])/freq_hi;
		if (fabs(pow_fft[j+1]) < 1e-30 && fabs(pow_fft[j]) < 1e-30)
			temp2[i]= 0.0;
		else
			temp2[i]= (180./M_PI)*atan2(pow_fft[j], pow_fft[j+1]);
	}
	for (i=0; i< freq_hi; i++) {
		pow_fft[i]= temp1[i];
		pow_fft[freq_hi+i]= temp2[i];
	} // }}}

	//*** Computing mag and phase of the voltage signal ** {{{

	temp1[0]= fabs(phi_fft[0])/freq_hi/2;
	temp2[0]= (phi_fft[0] > 0.0) ? 0.0 : 180.0;
	for (i=1, j=2; i< freq_hi; i++, j+=2) {
		temp1[i]= sqrt(phi_fft[j]*phi_fft[j]+ phi_fft[j+1]*phi_fft[j+1])/freq_hi;
		if (fabs(phi_fft[j+1]) < 1e-30 && fabs(phi_fft[j]) < 1e-30)
			temp2[i]= 0.0;
		else
			temp2[i]= (180./M_PI)*atan2(phi_fft[j], phi_fft[j+1]);
	}
	for (i=0; i< freq_hi; i++) {
		phi_fft[i]= temp1[i];
		phi_fft[freq_hi+i]= temp2[i];
	} // }}}

	//*** Computing mag and phase of the mid-potential signal ** {{{

	temp1[0]= fabs(mphi_fft[0])/freq_hi/2;
	temp2[0]= (mphi_fft[0] > 0.0) ? 0.0 : 180.0;
	for (i=1, j=2; i< freq_hi; i++, j+=2) {
		temp1[i]= sqrt(mphi_fft[j]*mphi_fft[j]+ mphi_fft[j+1]*mphi_fft[j+1])/freq_hi;
		if (fabs(mphi_fft[j+1]) < 1e-30 && fabs(mphi_fft[j]) < 1e-30)
			temp2[i]= 0.0;
		else
			temp2[i]= (180./M_PI)*atan2(mphi_fft[j], mphi_fft[j+1]);
	}
	for (i=0; i< freq_hi; i++) {
		mphi_fft[i]= temp1[i];
		mphi_fft[freq_hi+i]= temp2[i];
	} // }}}
} // }}}


// Modeline for ViM {{{1
// vim:set ts=4:
// vim600:fdm=marker fdl=0 fdc=3:
// 1}}}