C / C++ / MFC

<br />
void ei(double *x, double *expi, int nSize)<br />
{<br />
	enum{nSize = 100};<br />
	const int MAXIT = 100;<br />
	const double EULER=0.577215664901533;<br />
	const double EPS=1e-20;		//note arbitrary value of epsilon<br />
<br />
/*MAXIT is the maximum number of iterations allowed (for speed); <br />
	EULER is Euler's constant, EPS is epsilon, the relative or absolute<br />
	error; */<br />
<br />
	int k;<br />
	double fact[nSize], prev[nSize],sum[nSize],term[nSize];<br />
<br />
	for(int t=0;t<nSize;t++)<br />
	{<br />
		if (x[t]<=0.0)<br />
		{<br />
			AfxMessageBox("x value incorrect, in ei!", MB_OKCANCEL|MB_ICONSTOP);//THIS IS THE ERROR MSG I KEEP GETTING AT RUNTIME.<br />
			exit(1);<br />
		}<br />
<br />
		if (x[t]<EPS)<br />
		{<br />
			expi[t]=log(x[t])+EULER;	//Go ahead and converge!; since x<EPS, anyway<br />
		}<br />
<br />
		if (x[t]<=-log(EPS))<br />
		{<br />
			sum[0]=0.0;<br />
			fact[0]=1.0;<br />
	<br />
			for (k=1;k<MAXIT;k++)<br />
			{<br />
				fact[t] *= x[t]/k;<br />
				term[t]=fact[t]/k;<br />
				sum[t] += term[t];<br />
				if (term[t]<EPS*sum[t]) break;<br />
			}<br />
	<br />
			if (k>MAXIT)<br />
			{<br />
				AfxMessageBox("Series failed in ei!",MB_OKCANCEL|MB_ICONSTOP);<br />
				exit(2);<br />
			}<br />
	<br />
			expi[t]=sum[t]+log(x[t])+EULER;<br />
		}<br />
<br />
		else<br />
		{<br />
			sum[0]=0.0;<br />
			term[0]=1.0;<br />
<br />
			for (k=1;k<MAXIT;k++)<br />
			{<br />
				prev[t]=term[t];<br />
				term[t] *= k/x[t];<br />
<br />
				if (term[t]<EPS) break;<br />
				//since final sum is greater than 1, term itself approximates the relative error<br />
<br />
				if (term[t]<prev[t]) sum[t] += term[t];<br />
<br />
				else<br />
				{<br />
					sum[t] -= prev[t];<br />
					break;<br />
				}<br />
			}<br />
			<br />
			expi[t]=exp(x[t])*(1.0*sum[t])/x[t];<br />
		}<br />
	}<br />
}

<br />
void unbound(double *unbd, int nSize)<br />
{<br />
	enum{nSize = 100};<br />
	double A[nSize],B[nSize],k0[nSize],k1[nSize],s[nSize],expi[nSize],a[nSize],b[nSize];<br />
	<br />
	//trial: K0(x) --> x = rw*sqrt(s), as per A-32<br />
	s[0] = 1;<br />
	for (int t=1;t<nSize;t++)<br />
	{<br />
		double x;//this variable has been initialized in my dialog-control as 1; it is changed by the user at run-time, anyway, and this change is captured in the next 2 lines of code, in the instance of the class, below...<br />
<br />
		CKuchukDlg MyClas;	<br />
		MyClas.m_dWellRad = x;<br />
<br />
		void bessk0(double *x,double *k0,int nSize);<br />
		void bessk1(double *x,double *k1, int nSize);<br />
		void ei(double *x, double *expi, int nSize);<br />
<br />
		s[t] = 2*t*PI;<br />
		a[t] = x*sqrt(s[t]);<br />
		b[t] = 2*sqrt(s[t]);<br />
<br />
		bessk0(a,k0,nSize);<br />
		bessk1(a,k1,nSize);<br />
		ei(b,expi,nSize);<br />
<br />
		A[t] = k0[t]/(x*sqrt(s[t])*k1[t]);<br />
		B[t] = (1-exp(-2*sqrt(s[t])))/(2*sqrt(s[t]));<br />
<br />
		unbd[t]=0.5*(A[t]-B[t]-expi[t]);<br />
	}<br />
}<br />