/* File SurfaceParam.c, written by Steve Andrews, 2008.
This code is in the public domain.  It is not copyrighted and may not be copyrighted. */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "SurfaceParam.h"

/* Comment out either or both of the following include lines if the library
files are unavailable */
#include "math2.h"
#include "random2.h"


/* Declarations for functions that are only used internally */
double interpolate1D(double *xdata,double *ydata,int n,double x);
double interpolate2D(double *xdata,double *ydata,double *zdata,int nx,int ny,double x,double y);


/******************************************************************************/
/************  THINGS FROM math2 AND random2 LIBRARY FILES  *******************/
/******************************************************************************/

#ifndef __math2_h
	#include <float.h>
	#define SQRT2 1.41421356237
	#define SQRTPI 1.7724538509
	#define SQRT2PI 2.50662827462

	double gammalnD(double x);
	double gammpD(double a,double x);
	double erfccD(double x);
	double experfcD(double x);
	double erfnD(double x);
	int iseven(int x);
	void linefitD(double *x,double *y,int n,double *m,double *b);
#endif


#ifndef __random2_h
	double randCCD(void);
	double randCCD(void) {
		return (double)rand()/RAND_MAX; }
#endif


/******************************************************************************/
/************************  FUNCTIONS FOR EXTERNAL USE  ************************/
/******************************************************************************/


/* surfaceprob. */
double surfaceprob(double k1,double k2,double dt,double difc,double *p2ptr,enum SurfParamAlgo algo) {
	double step,p1,p2,kapp,kp,p1lo,p1hi,c1,c2,caideal,casim,kapfp,kapbp;
	int it;
	const int itmax=16;

	step=sqrt(2.0*difc*dt);
	p1=p2=0;

	if(algo==SPAirrTrans) {																				// IrrTrans
		return surfaceprob(k1,k2,dt,difc,p2ptr,SPAirrAds); }

	else if(algo==SPAirrTransT) {																	// IrrTransT
		return surfaceprob(k1,k2,dt,difc,p2ptr,SPAirrAdsT); }

	else if(algo==SPAirrTransQ) {																	// IrrTransQ
		return surfaceprob(k1,k2,dt,difc,p2ptr,SPAirrAdsQ); }

	else if(algo==SPArevTrans) {																	// RevTrans
		kapfp=k1*dt/step;
		kapbp=k2*dt/step;
		c1=kapfp+kapbp;
		c2=1.0/(c1*c1)*(2*c1-SQRTPI/SQRT2+SQRTPI/SQRT2*experfcD(SQRT2*c1));
		p1=kapfp*c2;
		p2=kapbp*c2; }

	else if(algo==SPAirrAds || algo==SPAirrAdsT) {								// IrrAds, IrrAdsT
		kapp=k1*dt/step;
		p1=lookupirrevadsorb(kapp,1); }

	else if(algo==SPAirrAdsQ) {																		// IrrAdsQ
		kapp=k1*dt/step;
		if(kapp<=0) p1=0;
		else if(kapp>=0.9) p1=1.0;
		else {
			p1=kapp*(SQRT2PI+kapp*(-3.3332146+kapp*(3.356688+kapp*(-1.5209235))));
			if(p1>1.0) p1=1.0; }}

	else if(algo==SPAirrAdsEC) {																	// IrrAdsEC
		kapp=k1*dt/step;
		p1=kapp*SQRT2PI;
		if(p1<0) p1=0;
		else if(p1>1) p1=1; }

	else if(algo==SPArevAds) {																		// RevAds
		kapp=k1*dt/step;
		kp=k2*dt;
		p1=lookuprevads(kapp,kp,1,&p2); }

	else if(algo==SPArevAdsND) {																	// RevAdsND
		kapp=k1*dt/step;
		kp=k2*dt;
		p2=1.0-exp(-kp);
		caideal=kapp*step/kp;
		p1lo=p1=0;
		p1hi=1.0;
		for(it=0;it<itmax;it++) {
			p1=0.5*(p1lo+p1hi);
			casim=lookuprevadsorbnd(p1,p2);
			if(casim>caideal) p1hi=p1;
			else p1lo=p1; }
		p1=0.5*(p1lo+p1hi); }

	else if(algo==SPAirrDes) {																		// IrrDes
		if(k2<=0) k2=k1;
		p1=k1/k2*(1.0-exp(-k2*dt)); }

	else if(algo==SPArevDes) {																		// RevDes
		p2=surfaceprob(k2,k1,dt,difc,&p1,SPArevAds); }

	else if(algo==SPAirrFlip) {																		// IrrFlip
		if(k2<=0) k2=k1;
		p1=k1/k2*(1.0-exp(-k2*dt)); }

	else if(algo==SPArevFlip) {																		// RevFlip
		p1=k1/(k1+k2)*(1.0-exp(-(k1+k2)*dt));
		p2=k2/(k1+k2)*(1.0-exp(-(k1+k2)*dt)); }

	else {																												// unrecognized
		p1=p2=-1; }

	if(p2ptr) *p2ptr=p2;
	return p1; }


/* desorbdist */
double desorbdist(double step,enum SurfParamAlgo algo) {
	double x,y,y2;

	y=randCCD();
	y2=y*y;

	if(algo==SPAirrDes)											// irreversible desorption
		x=(0.571825*y-0.552246*y2)/(1.0-1.53908*y+0.546424*y2);
	else if(algo==SPAirrDesC)
		x=2.0/3.0*SQRT2/SQRTPI;
	else if(algo==SPArevAds)								// reversible adsorption and desorption
		x=(0.729614*y-0.70252*y2)/(1.0-1.47494*y+0.484371*y2);
	else if(algo==SPArevAdsC)
		x=1.0/2.0*SQRTPI/SQRT2;
	else																		// unrecognized algorithm
		x=-1;

	return step*x; }


/* surfacerate. */
double surfacerate(double p1,double p2,double dt,double difc,double *k2ptr,enum SurfParamAlgo algo) {
	double step,k1,k2,k1p,k2p,casim;

	step=sqrt(2.0*difc*dt);
	k1=k2=0;

	if(algo==SPAirrTrans)																					// IrrTrans
		return surfacerate(p1,p2,dt,difc,k2ptr,SPAirrAds);

	else if(algo==SPAirrTransT)																		// IrrTransT
		return surfacerate(p1,p2,dt,difc,k2ptr,SPAirrAdsT);

	else if(algo==SPAirrTransQ)																		// IrrTransQ
		return surfacerate(p1,p2,dt,difc,k2ptr,SPAirrAdsQ);

	else if(algo==SPArevTrans) {																	// RevTrans
		k1p=lookuprevtrans(p1,p2,&k2);
		k1=k1p*step/dt; }

	else if(algo==SPAirrAds || algo==SPAirrAdsT) {								// IrrAds, IrrAdsT
		k1p=lookupirrevadsorb(p1,0);
		k1=k1p*step/dt; }

	else if(algo==SPAirrAdsQ) {																		// IrrAdsQ
		if(p1<=0) k1p=0;
		else if(p1>=1) k1p=0.85797628;
		else k1p=p1*(1.0/SQRT2PI+p1*(0.24761325+p1*(0.00616431+p1*(0.20383929))));
		k1=k1p*step/dt; }

	else if(algo==SPAirrAdsEC)	{																	// IrrAdsEC
		if(p1<0) p1=0;
		else if(p1>1) p1=1;
		k1p=p1/SQRT2PI;
		k1=k1p*step/dt; }

	else if(algo==SPArevAds) {																		// RevAds
		k1p=lookuprevads(p1,p2,0,&k2p);
		k1=k1p*step/dt;
		k2=k2p/dt; }

	else if(algo==SPArevAdsND) {																	// RevAdsND
		k2p=-log(1.0-p2);
		casim=lookuprevadsorbnd(p1,p2);
		k1p=casim*k2p/step;
		k1=k1p*step/dt;
		k2=k2p/dt; }

	else if(algo==SPAirrDes) {																		// IrrDes
		if(p2==0) p2=p1;
		k1p=-p1/p2*log(1.0-p2);
		k1=k1p/dt; }

	else if(algo==SPArevDes) {																		// RevDes
		k2=surfacerate(p2,p1,dt,difc,&k1,SPArevAds); }

	else if(algo==SPAirrFlip) {																		// IrrFlip
		if(p2==0) p2=p1;
		k1p=-p1/p2*log(1.0-p2);
		k1=k1p/dt; }

	else if(algo==SPArevFlip) {																		// RevFlip
		k1p=-p1/(p1+p2)*log(1.0-p1-p2);
		k2p=-p2/(p1+p2)*log(1.0-p1-p2);
		k1=k1p/dt;
		k2=k2p/dt; }

	else {																												// unrecognized algorithm
		k1=k2=-1; }

	if(k2ptr) *k2ptr=k2;
	return k1; }



/******************************************************************************/
/**********************  PARAMETER CALCULATION FUNCTIONS  *********************/
/******************************************************************************/


/* interpolate1D. */
double interpolate1D(double *xdata,double *ydata,int n,double x) {
	int i;
	double *xd,ans;

	if(n<4) return -1;
	for(i=0;i<n && xdata[i]<x;i++);
	if(i<2) i=2;
	else if(i>n-2) i=n-2;

	xd=xdata+(i-2);
	ans=ydata[i-2]* (x-xd[1])*(x-xd[2])*(x-xd[3])/((xd[0]-xd[1])*(xd[0]-xd[2])*(xd[0]-xd[3]));
	ans+=ydata[i-1]*(x-xd[0])*(x-xd[2])*(x-xd[3])/((xd[1]-xd[0])*(xd[1]-xd[2])*(xd[1]-xd[3]));
	ans+=ydata[i]*  (x-xd[0])*(x-xd[1])*(x-xd[3])/((xd[2]-xd[0])*(xd[2]-xd[1])*(xd[2]-xd[3]));
	ans+=ydata[i+1]*(x-xd[0])*(x-xd[1])*(x-xd[2])/((xd[3]-xd[0])*(xd[3]-xd[1])*(xd[3]-xd[2]));
	return ans; }


/* y is fast-changing index of zdata */
double interpolate2D(double *xdata,double *ydata,double *zdata,int nx,int ny,double x,double y) {
	int i,j;
	double ans,*xd,*yd,xcoeff[4],zval[4];

	if(nx<4 || ny<4) return -1;
	for(i=0;i<nx && xdata[i]<x;i++);		// i is index for x
	for(j=0;j<ny && ydata[j]<y;j++);		// j is index for y

	if(i<2) i=2;
	else if(i>nx-2) i=nx-2;
	if(j<2) j=2;
	else if(j>ny-2) j=ny-2;

	xd=xdata+(i-2);
	xcoeff[0]=(x-xd[1])*(x-xd[2])*(x-xd[3])/((xd[0]-xd[1])*(xd[0]-xd[2])*(xd[0]-xd[3]));
	xcoeff[1]=(x-xd[0])*(x-xd[2])*(x-xd[3])/((xd[1]-xd[0])*(xd[1]-xd[2])*(xd[1]-xd[3]));
	xcoeff[2]=(x-xd[0])*(x-xd[1])*(x-xd[3])/((xd[2]-xd[0])*(xd[2]-xd[1])*(xd[2]-xd[3]));
	xcoeff[3]=(x-xd[0])*(x-xd[1])*(x-xd[2])/((xd[3]-xd[0])*(xd[3]-xd[1])*(xd[3]-xd[2]));

	zval[0]=xcoeff[0]*zdata[ny*(i-2)+(j-2)]+xcoeff[1]*zdata[ny*(i-1)+(j-2)]+xcoeff[2]*zdata[ny*(i)+(j-2)]+xcoeff[3]*zdata[ny*(i+1)+(j-2)];
	zval[1]=xcoeff[0]*zdata[ny*(i-2)+(j-1)]+xcoeff[1]*zdata[ny*(i-1)+(j-1)]+xcoeff[2]*zdata[ny*(i)+(j-1)]+xcoeff[3]*zdata[ny*(i+1)+(j-1)];
	zval[2]=xcoeff[0]*zdata[ny*(i-2)+(j-0)]+xcoeff[1]*zdata[ny*(i-1)+(j-0)]+xcoeff[2]*zdata[ny*(i)+(j-0)]+xcoeff[3]*zdata[ny*(i+1)+(j-0)];
	zval[3]=xcoeff[0]*zdata[ny*(i-2)+(j+1)]+xcoeff[1]*zdata[ny*(i-1)+(j+1)]+xcoeff[2]*zdata[ny*(i)+(j+1)]+xcoeff[3]*zdata[ny*(i+1)+(j+1)];

	yd=ydata+(j-2);
	ans=zval[0]* (y-yd[1])*(y-yd[2])*(y-yd[3])/((yd[0]-yd[1])*(yd[0]-yd[2])*(yd[0]-yd[3]));
	ans+=zval[1]*(y-yd[0])*(y-yd[2])*(y-yd[3])/((yd[1]-yd[0])*(yd[1]-yd[2])*(yd[1]-yd[3]));
	ans+=zval[2]*(y-yd[0])*(y-yd[1])*(y-yd[3])/((yd[2]-yd[0])*(yd[2]-yd[1])*(yd[2]-yd[3]));
	ans+=zval[3]*(y-yd[0])*(y-yd[1])*(y-yd[2])/((yd[3]-yd[0])*(yd[3]-yd[1])*(yd[3]-yd[2]));

	return ans; }


/* lookupirrevadsorb */
double lookupirrevadsorb(double value,int pfromk) {
	/* data were generated with xdfmaketableirrev with n=401 and eps=0.0001. */
	double ponlist[]={
		0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1};
	double kappalist[]={
		0,0.0205318,0.0422852,0.0653621,0.089876,0.115953,0.143735,0.17338,0.205067,0.238998,0.275401,0.314537,0.356702,0.402241,0.451546,0.505076,0.563366,0.627042,0.696848,0.773667,0.858559};
	int npon=21;

	double ans;

	if(pfromk) {							// probability from kappa
		if(value<=0) ans=0;
		else if(value>=kappalist[npon-1]) ans=1.0;
		else ans=interpolate1D(kappalist,ponlist,npon,value); }
	else {										// kappa from probability
		if(value<=0) ans=0;
		else if(value>=1.0) ans=kappalist[npon-1];
		else ans=interpolate1D(ponlist,kappalist,npon,value); }
	return ans; }


/* lookuprevadsorbnd. */
double lookuprevadsorbnd(double probon,double proboff) {
	double ponlist[]={
		0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1};
	double pofflist[]= {
		0.02,0.07,0.12,0.17,0.22,0.27,0.32,0.37,0.42,0.47,0.52,0.57,0.62,0.67,0.72,0.77,0.82,0.87,0.92,0.97,1.02};
	/* poff is fast-changing index (columns), pon is slow-changing (rows). */
	/* data were generated with xdfmaketable with n=200 and eps=0.01. */
	double cstable[]={
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		1.01586,0.290586,0.169508,0.119653,0.0924592,0.0753371,0.0635657,0.0549757,0.048431,0.0432788,0.0391173,0.035686,0.0328081,0.0303597,0.0282514,0.0264169,0.0248061,0.0233805,0.0221098,0.0209701,0.020341,
		2.05784,0.592864,0.345837,0.24412,0.188638,0.153705,0.129689,0.112163,0.0988106,0.0882989,0.0798086,0.0728078,0.0669362,0.061941,0.0576395,0.0538967,0.0506103,0.0477017,0.0451092,0.042784,0.0415005,
		3.11512,0.904541,0.52956,0.373807,0.288851,0.23536,0.198585,0.171749,0.151303,0.135207,0.122206,0.111486,0.102496,0.0948466,0.08826,0.0825288,0.0774966,0.0730428,0.0690731,0.0655126,0.0635472,
		4.18328,1.22704,0.718315,0.509163,0.393444,0.320584,0.270493,0.23394,0.20609,0.184165,0.166457,0.151856,0.139609,0.129191,0.120219,0.112413,0.105558,0.0994917,0.0940845,0.0892348,0.0865577,
		5.25919,1.55934,0.914582,0.648033,0.502806,0.409693,0.345679,0.298965,0.263374,0.235356,0.212725,0.194065,0.178415,0.1651,0.153635,0.143659,0.134899,0.127146,0.120236,0.114038,0.110617,
		6.34277,1.90048,1.11799,0.791942,0.617367,0.50304,0.42444,0.367083,0.323383,0.288981,0.261194,0.238282,0.219066,0.202718,0.18864,0.176391,0.165635,0.156116,0.147631,0.140021,0.135821,
		7.42969,2.24994,1.32831,0.942077,0.730169,0.60102,0.50711,0.438582,0.38637,0.345267,0.312068,0.284693,0.261734,0.242202,0.225382,0.210747,0.197897,0.186523,0.176386,0.167294,0.162275,
		8.52325,2.60762,1.54578,1.09804,0.851479,0.695682,0.594064,0.513785,0.45262,0.404469,0.365578,0.33351,0.306614,0.283732,0.264028,0.246884,0.23183,0.218506,0.206631,0.19598,0.190101,
		9.61719,2.97182,1.76995,1.25982,0.977925,0.799258,0.676032,0.587051,0.522457,0.466876,0.421985,0.384968,0.353923,0.32751,0.304767,0.284977,0.2676,0.252221,0.238513,0.226219,0.219432,
		10.7216,3.34291,2.00122,1.42766,1.1097,0.907434,0.767729,0.66549,0.587876,0.532815,0.481583,0.439339,0.403908,0.373766,0.34781,0.325225,0.305394,0.287843,0.272199,0.258168,0.250423,
		11.8208,3.71906,2.23902,1.60139,1.24668,1.02021,0.863641,0.748673,0.66094,0.591924,0.537257,0.496927,0.456852,0.422759,0.3934,0.367855,0.345425,0.325573,0.307879,0.292009,0.283248,
		12.9267,4.10185,2.48305,1.78171,1.38898,1.13829,0.963851,0.836099,0.73801,0.660777,0.59829,0.546929,0.505625,0.474787,0.441816,0.413126,0.387936,0.365641,0.345769,0.327946,0.318107,
		14.042,4.48887,2.73296,1.96828,1.53729,1.261,1.06908,0.927461,0.819104,0.73356,0.663968,0.60671,0.558522,0.517781,0.483503,0.461336,0.433206,0.408309,0.386119,0.366215,0.355229,
		15.1488,4.88331,2.99035,2.16108,1.6914,1.38998,1.17915,1.02369,0.904305,0.810023,0.733472,0.670076,0.616706,0.57133,0.532291,0.498208,0.468704,0.443751,0.429215,0.40709,0.394878,
		16.2665,5.27705,3.25372,2.36011,1.85249,1.52399,1.29429,1.12458,0.994056,0.890541,0.806753,0.737013,0.678528,0.628335,0.585175,0.547682,0.514442,0.485057,0.458845,0.435273,0.419405,
		17.3854,5.67833,3.52276,2.56529,2.01876,1.66367,1.4147,1.23037,1.08836,0.975581,0.883838,0.807736,0.743912,0.689048,0.641601,0.600139,0.563809,0.531471,0.50218,0.475606,0.456787,
		18.4902,6.08248,3.7971,2.77652,2.19113,1.81047,1.54162,1.34214,1.18748,1.06513,0.965472,0.882717,0.812885,0.753151,0.701452,0.65624,0.616326,0.580766,0.548755,0.51939,0.497833,
		19.6254,6.48891,4.07638,2.99365,2.37171,1.96227,1.6735,1.45869,1.29256,1.16022,1.05162,0.961954,0.886212,0.821363,0.765188,0.715606,0.672193,0.633479,0.598584,0.566493,0.541924,
		20.7461,6.897,4.36015,3.22006,2.55677,2.12216,1.81286,1.58219,1.40235,1.25983,1.1434,1.04648,0.963866,0.893681,0.832827,0.779,0.731891,0.689372,0.651438,0.616458,0.588926,
		21.8668,7.31666,4.64797,3.44908,2.74795,2.28676,1.95716,1.71065,1.51915,1.36601,1.23969,1.13535,1.04699,0.970431,0.904691,0.846401,0.795412,0.749271,0.708115,0.669513,0.639327};
	const int npon=21,npoff=21;
	double zero=0;

	if(probon<=0) return 0;
	if(proboff<=0) return 1.0/zero;
	if(probon>1.0) probon=1.0;
	if(proboff>1.0) proboff=1.0;
	return interpolate2D(ponlist,pofflist,cstable,npon,npoff,probon,proboff); }


/* lookuprevads */
double lookuprevads(double value1,double value2,int pfromk,double *ans2ptr) {
	/* data were generated with Mathematica file RevAds.nb and were exported to file revadslookup.dat,
	kr is fast-changing index. */
	double xlist[]={
		0.01,0.059,0.108,0.157,0.206,0.255,0.304,0.353,0.402,0.451,0.5,0.549,0.598,0.647,0.696,0.745,0.794,0.843,0.892,0.941,0.99};
	double adsprobtable[]={
		0.0249247,0.0242857,0.0236019,0.0228691,0.022082,0.0212354,0.0203232,0.0193387,0.0182749,0.0171246,0.0158806,0.0145365,0.0130883,0.0115363,0.00988846,0.00816621,0.00641028,0.00468407,0.00306125,0.00158665,0.000255731,
		0.146554,0.142918,0.139023,0.134845,0.130352,0.125513,0.120292,0.114647,0.108537,0.101915,0.0947345,0.0869535,0.078539,0.0694807,0.0598105,0.0496336,0.039168,0.0287719,0.0188876,0.00982085,0.0015867,
		0.267179,0.260781,0.253922,0.246554,0.238624,0.230071,0.220827,0.210817,0.199959,0.188164,0.17534,0.161398,0.146263,0.129893,0.112315,0.0936803,0.0743426,0.05492,0.0362287,0.018905,0.00306256,
		0.386552,0.377647,0.36809,0.357814,0.346738,0.334774,0.321823,0.307773,0.2925,0.275868,0.257733,0.237949,0.216383,0.19294,0.167612,0.140556,0.112207,0.0833969,0.0553091,0.0289772,0.00470834,
		0.504415,0.493277,0.481311,0.468427,0.454522,0.439479,0.423166,0.405434,0.386113,0.365019,0.341947,0.316686,0.28903,0.258806,0.225939,0.190542,0.153074,0.114519,0.0763989,0.0402072,0.00655516,
		0.620501,0.607422,0.593355,0.578188,0.561794,0.544031,0.524732,0.503708,0.480747,0.455607,0.42802,0.397698,0.364346,0.327694,0.287559,0.24396,0.197309,0.14866,0.0998263,0.0528055,0.00864224,
		0.734533,0.719825,0.703984,0.68688,0.668365,0.648267,0.626388,0.6025,0.576343,0.547618,0.515987,0.481076,0.442488,0.399828,0.352769,0.301178,0.245337,0.186269,0.125994,0.0670368,0.0110198,
		0.846233,0.830222,0.812955,0.794284,0.774037,0.752017,0.727996,0.701706,0.672838,0.641035,0.605885,0.56692,0.523624,0.475454,0.421901,0.362615,0.297655,0.227884,0.155399,0.0832382,0.0137528,
		0.955321,0.938351,0.920023,0.900172,0.878607,0.855106,0.829411,0.801216,0.770166,0.735842,0.697755,0.655336,0.607939,0.554848,0.495332,0.428758,0.354851,0.274157,0.188667,0.101845,0.0169276,
		1.06153,1.04395,1.02494,1.00432,0.981872,0.957356,0.930485,0.900919,0.868256,0.832018,0.791637,0.746439,0.695635,0.638316,0.573489,0.500169,0.417623,0.325888,0.226588,0.123432,0.0206607,
		1.16459,1.14678,1.12748,1.10651,1.08363,1.05859,1.03107,1.0007,0.967037,0.929545,0.887577,0.840354,0.786935,0.726202,0.656861,0.577502,0.486806,0.384063,0.270181,0.148768,0.0251133,
		1.26427,1.2466,1.22741,1.20652,1.18368,1.15863,1.13101,1.10044,1.06643,1.0284,0.985625,0.937216,0.882085,0.81889,0.746005,0.661527,0.563407,0.449918,0.320775,0.178907,0.0305158,
		1.36035,1.34319,1.32453,1.30417,1.28185,1.2573,1.23017,1.20003,1.16638,1.12858,1.08584,1.03718,0.981357,0.916812,0.841562,0.753148,0.648652,0.525012,0.380134,0.215338,0.0372082,
		1.45265,1.43639,1.41866,1.39927,1.37797,1.35447,1.32841,1.29937,1.2668,1.23005,1.18827,1.1404,1.08506,1.02046,0.944273,0.853442,0.744044,0.611344,0.450651,0.260216,0.0457147,
		1.54104,1.52603,1.50964,1.49167,1.47188,1.44998,1.42561,1.39835,1.36765,1.33282,1.293,1.24706,1.19352,1.13039,1.055,0.963697,0.851452,0.711517,0.535642,0.316789,0.0568882,
		1.62542,1.61203,1.59737,1.58126,1.56346,1.54371,1.52165,1.49688,1.46885,1.43689,1.40011,1.35737,1.30712,1.24723,1.17475,1.08547,0.973219,0.828972,0.63982,0.390163,0.0722144,
		1.70574,1.69431,1.68177,1.66794,1.65262,1.63557,1.61645,1.59489,1.57037,1.54225,1.50969,1.47155,1.42628,1.37171,1.30471,1.22064,1.11233,0.968346,0.770094,0.488801,0.0945326,
		1.78201,1.77286,1.76279,1.75167,1.7393,1.72548,1.70993,1.69231,1.67217,1.64894,1.62184,1.58984,1.55148,1.50466,1.44629,1.37155,1.27263,1.13603,0.93695,0.627703,0.130025,
		1.85429,1.84771,1.84046,1.83241,1.82345,1.81339,1.80203,1.78909,1.77422,1.75696,1.73669,1.71253,1.68325,1.64705,1.60115,1.54108,1.45916,1.34107,1.157,0.835808,0.195157,
		1.92266,1.91893,1.9148,1.9102,1.90506,1.89928,1.89271,1.8852,1.87652,1.86637,1.85436,1.83991,1.8222,1.79999,1.77131,1.73286,1.67866,1.59658,1.45791,1.17558,0.352706,
		1.98727,1.98661,1.98588,1.98507,1.98416,1.98313,1.98196,1.98061,1.97905,1.9772,1.975,1.97232,1.96901,1.96478,1.95921,1.95155,1.94034,1.92236,1.88884,1.80435,1.19212};
	double kappartable[]={
		0.00401057,0.00411327,0.00422333,0.00434171,0.00446954,0.00460816,0.00475924,0.00492481,0.00510743,0.00531035,0.0055378,0.0057954,0.00609085,0.00643501,0.00684395,0.00734275,0.00797365,0.00881587,0.0100456,0.012325,0.0199136,
		0.0236948,0.0243079,0.0249657,0.0256741,0.0264401,0.027272,0.0281804,0.0291779,0.0302809,0.03151,0.0328926,0.0344652,0.0362785,0.0384052,0.0409551,0.0441045,0.0481614,0.0537388,0.0623447,0.0802753,0.0944274,
		0.0434378,0.0445739,0.0457943,0.0471102,0.0485351,0.0500853,0.0517811,0.0536477,0.0557169,0.0580301,0.0606418,0.0636264,0.0670879,0.0711789,0.0761343,0.0823432,0.0905177,0.102177,0.121593,0.157299,0.131216,
		0.0632464,0.0649191,0.0667181,0.0686605,0.0707669,0.0730628,0.0755793,0.0783558,0.0814424,0.0849045,0.0888296,0.0933379,0.0986006,0.104873,0.11256,0.122355,0.135595,0.155395,0.19227,0.236135,0.160189,
		0.0831276,0.0853517,0.0877466,0.0903361,0.0931489,0.0962203,0.0995941,0.103326,0.107487,0.112172,0.117506,0.123667,0.130912,0.139629,0.150456,0.164527,0.184173,0.215466,0.285302,0.250177,0.184053,
		0.103089,0.10588,0.10889,0.112149,0.115695,0.119575,0.123847,0.128584,0.133884,0.139874,0.146728,0.154693,0.164134,0.175618,0.190102,0.209376,0.237409,0.286141,0.43914,0.259578,0.204436,
		0.123139,0.126515,0.130159,0.134112,0.138421,0.143145,0.14836,0.154159,0.16067,0.16806,0.176561,0.186507,0.198402,0.213051,0.231861,0.25762,0.297106,0.375162,0.348514,0.271392,0.222314,
		0.143286,0.147264,0.151566,0.15624,0.161344,0.166952,0.173158,0.180083,0.187886,0.196784,0.207079,0.219216,0.233877,0.25219,0.27621,0.310283,0.366288,0.501817,0.35197,0.283444,0.238317,
		0.163538,0.16814,0.173123,0.178546,0.18448,0.191017,0.19827,0.206391,0.215578,0.226108,0.23837,0.252947,0.270757,0.293371,0.323792,0.368896,0.450479,0.727019,0.354706,0.295144,0.25287,
		0.183906,0.189153,0.194843,0.201047,0.207851,0.215364,0.223727,0.233122,0.243799,0.256104,0.270534,0.287851,0.309286,0.337025,0.375501,0.435864,0.56093,0.430173,0.35966,0.306346,0.266273,
		0.2044,0.210314,0.216741,0.223761,0.231477,0.240021,0.249561,0.260321,0.272608,0.286853,0.303692,0.324113,0.349769,0.383724,0.432619,0.515193,0.723369,0.430663,0.365873,0.317042,0.278747,
		0.22503,0.231639,0.238831,0.246705,0.255382,0.265016,0.275811,0.288038,0.302072,0.318453,0.337985,0.361956,0.392594,0.434239,0.497053,0.614007,1.00754,0.431193,0.372718,0.327267,0.290456,
		0.245809,0.253139,0.261132,0.269902,0.279591,0.290382,0.302518,0.316327,0.332269,0.351012,0.373582,0.401658,0.438264,0.489636,0.57176,0.745819,0.497567,0.433218,0.379858,0.337067,0.301526,
		0.266748,0.27483,0.283661,0.293373,0.304132,0.316155,0.329731,0.345255,0.363289,0.384664,0.410687,0.443564,0.487441,0.551434,0.661511,0.937958,0.496691,0.436454,0.387118,0.346488,0.312055,
		0.287862,0.296728,0.306438,0.317143,0.329036,0.342375,0.357501,0.374892,0.395234,0.419564,0.449551,0.488115,0.541011,0.621852,0.774304,0.856578,0.496207,0.440523,0.394398,0.355573,0.322122,
		0.309164,0.318851,0.329484,0.341238,0.354338,0.369086,0.38589,0.405323,0.428226,0.455897,0.490481,0.535876,0.600188,0.704218,0.924072,0.555998,0.496815,0.445151,0.401645,0.364359,0.33179,
		0.330669,0.341218,0.352823,0.365687,0.380075,0.396339,0.414964,0.436643,0.462408,0.493888,0.533862,0.587585,0.666664,0.80365,1.10195,0.554415,0.498436,0.450153,0.408829,0.372878,0.341108,
		0.352396,0.363848,0.37648,0.390525,0.40629,0.42419,0.444803,0.468965,0.497948,0.53381,0.580175,0.644219,0.74284,0.928211,1.02946,0.553377,0.500858,0.455406,0.415935,0.381159,0.35012,
		0.374362,0.386765,0.400483,0.415786,0.43303,0.452703,0.475495,0.502419,0.535048,0.575997,0.630037,0.707089,0.832175,1.08215,0.91603,0.553216,0.503888,0.460828,0.422955,0.389225,0.35886,
		0.396586,0.409992,0.424863,0.441509,0.460347,0.48195,0.507142,0.537156,0.573947,0.620862,0.684241,0.777971,0.939742,1.1105,0.922459,0.553899,0.507377,0.466361,0.429885,0.397098,0.367357,
		0.419089,0.433557,0.449655,0.467741,0.4883,0.512012,0.539862,0.573356,0.614938,0.66892,0.743813,0.859305,1.09787,1.05978,0.604705,0.555297,0.511213,0.471966,0.436726,0.404794,0.375635};
	double krtable[]={
		0.00999264,0.0575665,0.102992,0.146473,0.188198,0.228343,0.267075,0.304553,0.340934,0.376377,0.411046,0.445115,0.478784,0.512286,0.545915,0.580071,0.615344,0.65271,0.694086,0.929393,0.920201,
		0.0102059,0.0587783,0.105134,0.149489,0.192041,0.232977,0.272472,0.310697,0.34782,0.38401,0.419445,0.45432,0.488855,0.52332,0.55806,0.593559,0.630566,0.670398,0.71589,0.929136,0.872188,
		0.0104296,0.0600492,0.107381,0.152651,0.196071,0.237837,0.278136,0.31715,0.355057,0.392042,0.428298,0.464042,0.499523,0.535051,0.571038,0.608084,0.647156,0.690084,0.74132,0.926224,0.864545,
		0.0106646,0.0613838,0.109739,0.155972,0.200303,0.242943,0.28409,0.323937,0.362677,0.40051,0.437649,0.474334,0.510851,0.547562,0.584965,0.623812,0.665389,0.712337,0.77221,0.912064,0.86288,
		0.0109117,0.0627872,0.112219,0.159463,0.204754,0.248315,0.290357,0.331087,0.370714,0.409455,0.447546,0.485257,0.522917,0.560955,0.59998,0.640957,0.685642,0.738028,0.812481,0.887084,0.86304,
		0.0111721,0.0642651,0.11483,0.163139,0.209442,0.253975,0.296964,0.338633,0.379206,0.418923,0.458045,0.49688,0.535809,0.575347,0.616254,0.659793,0.70844,0.768603,0.87286,0.879146,0.863976,
		0.0114467,0.0658238,0.117584,0.167017,0.214388,0.259949,0.303944,0.346612,0.388199,0.428968,0.469213,0.509285,0.549633,0.590885,0.634003,0.680685,0.734559,0.806705,0.911966,0.876102,0.865256,
		0.0117369,0.0674704,0.120493,0.171113,0.219614,0.266266,0.31133,0.355066,0.397742,0.43965,0.481123,0.522567,0.564517,0.607747,0.653505,0.704134,0.765195,0.857747,0.897971,0.874887,0.866658,
		0.0120441,0.0692128,0.123571,0.175448,0.225147,0.272958,0.319162,0.364042,0.407893,0.451041,0.493864,0.53684,0.580614,0.626157,0.675121,0.730846,0.802317,0.93333,0.890728,0.874495,0.868052,
		0.0123699,0.0710601,0.126835,0.180044,0.231017,0.280062,0.327485,0.373595,0.418718,0.463221,0.507541,0.55224,0.598111,0.646397,0.699335,0.761857,0.849347,0.915924,0.886683,0.874472,0.869352,
		0.0127159,0.0730224,0.130301,0.184928,0.237256,0.28762,0.336351,0.383788,0.430295,0.476288,0.522276,0.568932,0.617244,0.668832,0.726812,0.798751,0.912348,0.905637,0.884117,0.874574,0.870495,
		0.0130844,0.0751112,0.133991,0.190129,0.243904,0.29568,0.345818,0.394693,0.442714,0.490355,0.538216,0.587116,0.638304,0.693941,0.758487,0.844015,1.00208,0.89859,0.882255,0.874649,0.871436,
		0.0134776,0.0773395,0.137927,0.195679,0.251003,0.304297,0.355956,0.406396,0.456079,0.505556,0.555538,0.607039,0.661666,0.722367,0.795717,0.901603,0.918532,0.893477,0.880697,0.874599,0.872137,
		0.0138981,0.0797225,0.142138,0.201617,0.258605,0.313535,0.366843,0.418994,0.470516,0.522051,0.574461,0.629012,0.687817,0.75499,0.840503,0.977596,0.909149,0.889378,0.879206,0.874349,0.872568,
		0.0143491,0.0822773,0.146652,0.207987,0.266768,0.323468,0.378573,0.432605,0.486173,0.540037,0.59525,0.653427,0.717398,0.793041,0.895829,1.02959,0.90143,0.885763,0.877629,0.873847,0.872703,
		0.0148339,0.085024,0.151507,0.214842,0.27556,0.334184,0.391255,0.447368,0.503228,0.559751,0.61824,0.68079,0.751269,0.838264,0.966102,0.916444,0.894791,0.88231,0.875861,0.873047,0.872516,
		0.0153568,0.0879858,0.156743,0.22224,0.28506,0.345784,0.405019,0.463447,0.5219,0.581486,0.643851,0.711759,0.790606,0.893163,1.05622,0.906282,0.888722,0.878811,0.873826,0.871914,0.871989,
		0.0159226,0.0911899,0.162409,0.230252,0.295363,0.35839,0.42002,0.481044,0.542452,0.605608,0.672623,0.747212,0.837037,0.961352,1.06485,0.896898,0.882869,0.875128,0.871463,0.870417,0.8711,
		0.0165368,0.0946684,0.168563,0.238962,0.306581,0.372147,0.436446,0.500404,0.565213,0.63258,0.705257,0.788325,0.89285,1.04702,1.04351,0.887988,0.876988,0.871157,0.868726,0.868529,0.869832,
		0.0172062,0.0984591,0.175272,0.248469,0.318847,0.387229,0.454524,0.521826,0.590594,0.662995,0.742678,0.836714,0.961304,1.07521,1.04595,0.879253,0.870913,0.866822,0.865574,0.866227,0.868168,
		0.0179387,0.102608,0.182619,0.258892,0.332323,0.403851,0.474535,0.545687,0.619117,0.697623,0.786125,0.894623,1.06788,1.08185,0.881251,0.870459,0.864519,0.862065,0.861976,0.863488,0.866092};
	const int nx=21;
	double ans1,ans2,kappar,kr;

	if(pfromk) {
		if(value1<0) value1=0;
		if(value2<0) value2=0;
		kappar=value1/(1.0+value1);
		kr=value2/(1.0+value2);
		ans1=interpolate2D(xlist,xlist,adsprobtable,nx,nx,kappar,kr);
		if(ans1<0) ans1=0;
		else if(ans1>1) ans1=1;
		ans2=ans1*value2/value1/SQRT2PI; }
	else {
		if(value1<0) value1=0;
		else if(value1>1) value1=1;
		if(value2<0) value2=0;
		else if(value2>1) value2=1;
		kappar=interpolate2D(xlist,xlist,kappartable,nx,nx,value1,value2);
		kr=interpolate2D(xlist,xlist,krtable,nx,nx,value1,value2);
		if(kappar<0) kappar=0;
		else if(kappar>0.999999999) kappar=0.999999999;
		if(kr<0) kr=0;
		else if(kr>0.999999999) kr=0.999999999;
		ans1=kappar/(1.0-kappar);
		ans2=kr/(1.0-kr); }

	if(ans2ptr) *ans2ptr=ans2;
	return ans1; }



/* lookuprevtrans */
double lookuprevtrans(double pf,double pb,double *kbptr) {
	/* data were generated with Mathematica file RevTrans.nb */
	double xlist[]={
		0.01,0.059,0.108,0.157,0.206,0.255,0.304,0.353,0.402,0.451,0.5,0.549,0.598,0.647,0.696,0.745,0.794,0.843,0.892,0.941,0.99};
	double kappatable[]={
		0.00402361,0.0041102,0.00420107,0.00429655,0.00439701,0.00450288,0.0046146,0.0047327,0.00485776,0.00499044,0.00513147,0.0052817,0.00544209,0.00561374,0.00579792,0.00599608,0.00620994,0.00644149,0.00669308,0.00696751,0.00726812,
		0.0242502,0.0247863,0.0253497,0.0259424,0.026567,0.0272261,0.0279229,0.0286608,0.0294436,0.0302757,0.031162,0.0321083,0.0331211,0.0342077,0.0353769,0.0366386,0.0380048,0.0394892,0.0411083,0.0428819,0.0448337,
		0.0453716,0.0464028,0.0474878,0.0486311,0.0498376,0.0511131,0.0524638,0.0538967,0.0554199,0.0570424,0.0587746,0.0606284,0.0626175,0.0647577,0.0670673,0.0695681,0.0722853,0.0752491,0.0784957,0.0820684,0.0860203,
		0.0674559,0.0690331,0.0706951,0.0724492,0.0743034,0.0762668,0.0783498,0.0805641,0.0829227,0.0854409,0.0881358,0.0910273,0.0941385,0.0974961,0.101131,0.105081,0.10939,0.114109,0.119303,0.125048,0.131439,
		0.0905785,0.0927592,0.0950607,0.0974936,0.10007,0.102803,0.105708,0.108803,0.112107,0.115643,0.119437,0.123519,0.127925,0.132695,0.137878,0.143531,0.149723,0.156538,0.164076,0.172461,0.181848,
		0.114823,0.117672,0.120684,0.123873,0.127256,0.130852,0.134683,0.138773,0.14315,0.147847,0.1529,0.158353,0.164258,0.170674,0.177672,0.185337,0.193773,0.203103,0.213483,0.225103,0.238205,
		0.140284,0.143874,0.147676,0.151709,0.155997,0.160564,0.16544,0.170658,0.176257,0.182281,0.188782,0.195821,0.20347,0.211812,0.220951,0.231007,0.242131,0.254505,0.268359,0.283978,0.30173,
		0.167064,0.171479,0.176162,0.181141,0.186444,0.192106,0.198165,0.204666,0.211662,0.219211,0.227385,0.236266,0.245953,0.256565,0.268242,0.281159,0.295528,0.311614,0.329751,0.350364,0.374009,
		0.195282,0.200616,0.206285,0.212324,0.218772,0.225672,0.233076,0.241042,0.24964,0.258948,0.269062,0.280094,0.292178,0.305477,0.320187,0.33655,0.354869,0.375523,0.398999,0.425926,0.457138,
		0.225069,0.231429,0.238205,0.245438,0.25318,0.261486,0.270423,0.280068,0.290511,0.301858,0.314235,0.327792,0.342711,0.359214,0.377572,0.398124,0.421296,0.447633,0.477842,0.512859,0.553949,
		0.256573,0.264085,0.272105,0.280687,0.289896,0.299804,0.310497,0.322074,0.334654,0.348376,0.363406,0.379946,0.398242,0.418595,0.441379,0.467069,0.496267,0.529758,0.56858,0.614134,0.66836,
		0.289965,0.298771,0.308194,0.318306,0.329185,0.340926,0.353638,0.36745,0.382516,0.39902,0.417181,0.43727,0.459617,0.484634,0.512841,0.544901,0.581674,0.624301,0.674319,0.733859,0.80596,
		0.325437,0.335702,0.346715,0.358566,0.371354,0.385201,0.400246,0.416657,0.434633,0.454416,0.476297,0.500639,0.52789,0.558614,0.593535,0.633591,0.680022,0.734504,0.799359,0.877894,0.974993,
		0.363209,0.375125,0.387946,0.401783,0.416764,0.433042,0.450798,0.470247,0.49165,0.515325,0.541661,0.571145,0.604387,0.642169,0.685507,0.735742,0.794689,0.864858,0.949829,1.05488,1.18816,
		0.403535,0.417327,0.432212,0.448328,0.465839,0.484939,0.505861,0.528885,0.554353,0.582684,0.6144,0.65016,0.690803,0.737423,0.791463,0.854874,0.930357,1.02176,1.13478,1.27815,1.46609,
		0.446708,0.462641,0.479891,0.498635,0.51908,0.541475,0.56612,0.59338,0.623706,0.657655,0.695932,0.739438,0.78934,0.847184,0.915059,0.995857,1.0937,1.21467,1.36813,1.56931,1.84469,
		0.493069,0.511454,0.531431,0.55322,0.577088,0.603355,0.632408,0.664728,0.70091,0.741705,0.788072,0.841256,0.902905,0.975244,1.06136,1.16563,1.29456,1.45811,1.67253,1.96601,2.39241,
		0.543018,0.564227,0.587361,0.612702,0.640589,0.671436,0.70575,0.744165,0.787478,0.836706,0.893172,0.958626,1.03543,1.12686,1.23757,1.37445,1.5481,1.77574,2.08734,2.54005,3.25809,
		0.597023,0.621502,0.648316,0.677824,0.710464,0.746773,0.787421,0.833251,0.88534,0.945088,1.01435,1.09561,1.19235,1.3095,1.45434,1.63808,1.87896,2.20867,2.68769,3.44747,4.83783,
		0.655643,0.68393,0.715059,0.749491,0.787795,0.830676,0.879023,0.933974,0.997005,1.07007,1.1558,1.25785,1.38144,1.53423,1.72807,1.98218,2.33,2.83533,3.63685,5.10358,8.65281,
		0.719544,0.752293,0.788519,0.828817,0.873931,0.924796,0.982608,1.04892,1.12579,1.21598,1.32335,1.45337,1.61412,1.81805,2.08539,2.45133,2.98297,3.82623,5.36934,9.10338,31.1321};
	const int nx=21;
	double kf,kb;

	if(pf<0) pf=0;
	else if(pf>1) pf=1;
	if(pb<0) pb=0;
	else if(pb>1) pb=1;
	kf=interpolate2D(xlist,xlist,kappatable,nx,nx,pf,pb);
	if(kf<0) kf=0;

	if(kbptr) {
		kb=interpolate2D(xlist,xlist,kappatable,nx,nx,pb,pf);
		if(kb<0) kb=0;
		*kbptr=kb; }

	return kf; }



/******************************************************************************/
/********   FUNCTIONS FOR INVESTIGATING A PARTIALLY ADSORBING SURFACE  ********/
/******************************************************************************/


/* xdfdiffuse. */
void xdfdiffuse(double *x,double *xdfa,double *xdfd,int n) {
	int i,j;
	double x0,x1,f0,f1,sum,grn;

	for(i=0;i<n;i++) {
		x1=x[0];
		grn=1.0/SQRT2PI*exp(-(x[i]-x1)*(x[i]-x1)/2.0);
		f1=xdfa[0]*grn;
		sum=0;
		for(j=1;j<n;j++) {
			x0=x1;
			f0=f1;
			x1=x[j];
			grn=1.0/SQRT2PI*exp(-(x[i]-x1)*(x[i]-x1)/2.0);
			f1=xdfa[j]*grn;
			sum+=0.5*(f0+f1)*(x1-x0); }
		sum+=0.5*(1.0+erfnD((x[i]-x[n-1])/SQRT2));
		xdfd[i]=sum; }
	return; }


/* xdfadsorb. */
double xdfadsorb(double *x,double *xdf,int n,double probon) {
	int j,j2;
	double sum,x0,x1,f0,f1,cleft,reflect;

	x0=x1=x[0];						// measure total flux
	f0=f1=xdf[0];
	cleft=2.0*f1/(1.0+erfnD(x1/SQRT2));
	sum=cleft/2.0*(SQRT2/SQRTPI*exp(-x1*x1/2.0)+x1*(1.0+erfnD(x1/SQRT2)));
	for(j=1;x1<0;j++) {
		x0=x1;
		f0=f1;
		x1=x[j];
		f1=xdf[j];
		sum+=0.5*(f0+f1)*(x1-x0); }
	f0=0;
	sum-=0.5*(f0+f1)*(x1-x0);

	reflect=1.0-probon;		// perform reflection and adsorption
	if(reflect<0) reflect=0;
	j-=2;
	j2=j+1;
	for(;j>=0;j--) {
		xdf[j2++]+=reflect*xdf[j];
		xdf[j]=0; }

	return sum*probon; }


/* xdfdesorb. */
void xdfdesorb(double *x,double *xdf,int n,double b,double flux) {
	int i;
	double grn;

	for(i=0;i<n;i++) {
		grn=1.0/SQRT2PI*exp(-(x[i]-b)*(x[i]-b)/2.0);
		xdf[i]+=flux*grn; }
	return; }


/* xdfdesorbdelta. */
void xdfdesorbdelta(double *x,double *xdf,int n,double b,double flux) {
	int i;

	for(i=0;i<n-1 && x[i]<=b;i++);
	xdf[i]+=2.0*flux/(x[i+1]-x[i-1]);
	return; }


/* xdfsteadystate. */
double xdfsteadystate(double *x,double *xdfa,double *xdfd,int n,double cs,double b,double probon,double proboff,double eps) {
	const int maxit=100000;
	const double maxflux=1e7;
	int i,it;
	double fluxon,fluxoff,flux,fluxold;

	fluxon=flux=fluxold=0;
	fluxoff=cs*proboff;

	for(it=0;(it<30) || (it<floor(0.1/eps)) || (it<maxit && flux<maxflux && fabs((flux-fluxold)/(fluxold+1e-20))>eps);it++) {
		xdfdiffuse(x,xdfa,xdfd,n);
		if(proboff>0) xdfdesorb(x,xdfd,n,b,fluxoff);
		fluxon=xdfadsorb(x,xdfd,n,probon);
		fluxoff=cs*proboff;
		fluxold=flux;
		flux=fluxon-fluxoff;
		cs+=flux;
		for(i=0;i<n;i++) xdfa[i]=xdfd[i]; }

	xdfdesorbdelta(x,xdfd,n,b,fluxoff);
	if(it>=maxit || flux>=maxflux) cs=-1.0;
	if(proboff==0) return flux;
	return cs; }


/* xdfmaketableirrev. */
void xdfmaketableirrev(void) {
	const double ponlo=0,ponhi=1.0,ponincr=0.05;							// pon low, high, increment
	double *x,*xdfa,*xdfd,xlo,xhi,dx,flux1,flux2,slope1,slope2,intercept1,intercept2,probon,eps,xfitlo,xfithi;
	int i,i2,n,ifitlo,ifithi,npon;
	char ynmro[256],ynxdf[256];

	fprintf(stderr,"Enter the number of position points for the concentration (e.g. 200): ");
	scanf("%i",&n);
	if(iseven(n)) n++;
	fprintf(stderr,"Enter low and high x values (e.g. -6 and 10): ");
	scanf("%lf %lf",&xlo,&xhi);
	fprintf(stderr,"Enter fit domain for x values (e.g. 3 and 7): ");
	scanf("%lf %lf",&xfitlo,&xfithi);
	fprintf(stderr,"Enter epsilon (e.g. 0.0001): ");
	scanf("%lf",&eps);
	fprintf(stderr,"Do you want machine readable output (y/n)? ");
	scanf("%s",ynmro);
	if(ynmro[0]!='y') {
		fprintf(stderr,"Do you want xdf output (y/n)? ");
		scanf("%s",ynxdf); }
	else ynxdf[0]='n';

	x=(double*)calloc(n,sizeof(double));
	xdfa=(double*)calloc(n,sizeof(double));
	xdfd=(double*)calloc(n,sizeof(double));
	if(!x || !xdfa || !xdfd) {
		fprintf(stderr,"Out of memory.  Function stopped.\n");return; }

	dx=(xhi-xlo)/n;
	x[0]=xlo;
	for(i=1;x[i-1]<0 && i<n;i++) x[i]=x[i-1]+dx;
	x[i-1]=-0.0001;																// x point below 0 is set to -0.0001
	i2=i-1;
	for(;i2>=0 && i<n;i++) x[i]=-x[i2--];
	for(;i<n;i++) x[i]=x[i-1]+dx;

	ifitlo=ifithi=-1;
	for(i=0;i<n && ifitlo==-1;i++) if(x[i]>=xfitlo) ifitlo=i;
	for(;i<n && ifithi==-1;i++) if(x[i]>xfithi) ifithi=i-1;
	if(ifitlo==-1 || ifithi==-1) {
		fprintf(stderr,"Fit domain is not within x range\n");return; }

	if(ynmro[0]=='y') {
		npon=0;
		printf("\tconst double ponlist[]={\n\t\t");
		for(probon=ponlo;probon<ponhi+ponincr/2.0;probon+=ponincr) {
			npon++;
			printf("%g,",probon); }
		printf("};\n");
		printf("\tconst int npon=%i;\n",npon);
		printf("\t/* data were generated with xdfmaketableirrev with n=%i and eps=%g. */\n",n,eps);
		printf("\tconst double irrevtable[]={\n\t\t"); }
	else
		printf("P_ads flux1 flux2 slope1 slope2 inter1 inter2 K'\n");

	for(probon=ponlo;probon<ponhi+ponincr/2;probon+=ponincr) {
		for(i=0;i<n;i++) xdfa[i]=xdfd[i]=(x[i]<0?0:1);
		flux1=xdfsteadystate(x,xdfa,xdfd,n,0,0,probon,0,eps);
		linefitD(x+ifitlo,xdfa+ifitlo,ifithi-ifitlo+1,&slope1,&intercept1);
		for(i=0;i<n;i++) xdfa[i]=xdfd[i]=0;
		flux2=xdfsteadystate(x,xdfa,xdfd,n,0,0,probon,0,eps);
		linefitD(x+ifitlo,xdfa+ifitlo,ifithi-ifitlo+1,&slope2,&intercept2);
		if(ynmro[0]=='y') printf("%g,",(flux1+flux2)/(intercept1+intercept2));
		else printf("%g %g %g %g %g %g %g %g\n",probon,flux1,flux2,slope1,slope2,intercept1,intercept2,(flux1+flux2)/(intercept1+intercept2));
		if(ynxdf[0]=='y') {
			xdfdiffuse(x,xdfa,xdfd,n);
			for(i=0;i<n;i++) printf("%g %g %g\n",x[i],xdfa[i],xdfd[i]); }}
	if(ynmro[0]=='y') printf("};\n");

	return; }


/* xdfmaketable */
void xdfmaketable(void) {
	const double ponlo=0,ponhi=1.0,ponincr=0.05;							// pon low, high, increment
	const double pofflo=0.02,poffhi=1.0,poffincr=0.05;				// poff low, high, increment
	double eps,pon,kon,koff,poff,cs,*x,*xdfa,*xdfd,dx;
	int i,i2,n,npon,npoff;
	char yn[256];

	fprintf(stderr,"\nFunction for calculating steady-state surface concentrations\n");
	fprintf(stderr,"for various adsorption and desorption probabilities.\n\n");
	fprintf(stderr,"Enter the number of position points for the concentration (e.g. 200): ");
	scanf("%i",&n);
	if(n<10) {fprintf(stderr,"Value is too low.  Function stopped.\n");return; }
	if(iseven(n)) n++;
	fprintf(stderr,"Enter level of precision (e.g. 1e-4): ");
	scanf("%lf",&eps);
	if(eps<=0) {fprintf(stderr,"Impossible precision.  Function stopped.\n");return; }
	fprintf(stderr,"Do you want machine readable output (y/n)? ");
	scanf("%s",yn);

	x=(double*)calloc(n,sizeof(double));
	xdfa=(double*)calloc(n,sizeof(double));
	xdfd=(double*)calloc(n,sizeof(double));
	if(!x || !xdfa || !xdfd) {
		fprintf(stderr,"Out of memory.  Function stopped.\n");return; }

	dx=16.0/n;																		// x goes from -6 to about +10
	x[0]=-6.0;																		// x min is set to -6
	for(i=1;x[i-1]<0 && i<n;i++) x[i]=x[i-1]+dx;
	x[i-1]=-0.0001;																// x point below 0 is set to -0.0001
	i2=i-1;
	for(;i2>=0 && i<n;i++) x[i]=-x[i2--];
	for(;i<n;i++) x[i]=x[i-1]+dx;

	if(yn[0]=='y') {
		npon=npoff=0;
		printf("\tconst double ponlist[]={\n\t\t");
		for(pon=ponlo;pon<ponhi+ponincr/2.0;pon+=ponincr) {
			npon++;
			printf("%g,",pon); }
		printf("};\n");
		printf("\tconst double pofflist[]= {\n\t\t");
		for(poff=pofflo;poff<poffhi+poffincr/2.0;poff+=poffincr) {
			npoff++;
			printf("%g,",poff); }
		printf("};\n");
		printf("\tconst int npon=%i,npoff=%i;\n",npon,npoff);
		printf("\t/* poff is fast-changing index (columns), pon is slow-changing (rows). */\n");
		printf("\t/* data were generated with xdfmaketable with n=%i and eps=%g. */\n",n,eps);
		printf("\tconst double cstable[]={\n\t\t"); }

	for(pon=ponlo;pon<ponhi+ponincr/2.0;pon+=ponincr) {
		if(pon>ponhi) pon=ponhi;
		for(poff=pofflo;poff<poffhi+poffincr/2.0;poff+=poffincr) {
			if(poff>poffhi) poff=poffhi;
			for(i=0;i<n;i++) xdfa[i]=xdfd[i]=(x[i]<0?0:1);
			kon=pon/SQRT2PI;
			koff=-log(1.0-poff);
			cs=poff<1.0?kon/koff:0;
			cs=xdfsteadystate(x,xdfa,xdfd,n,cs,0.0,pon,poff,eps);
			if(yn[0]=='y') printf("%g,",cs);
			else printf("%g %g %g\n",pon,poff,cs); }
		if(yn[0]=='y') printf("\n\t\t"); }
	if(yn[0]=='y') printf("};\n");

	free(x);
	free(xdfa);
	free(xdfd);
	fprintf(stderr,"Done making table\n");

	return; }



/******************************************************************************/
/********************  CODE FROM math2 LIBRARY FILE  **************************/
/******************************************************************************/

#ifndef __math2_h

	double gammalnD(double x)	{	/* Returns natural log of gamma function, partly from Numerical Recipies and part from Math CRC */
		double sum,t;
		int j;
		static double c[6]={76.18009173,-86.50532033,24.01409822,-1.231739516,0.120858003e-2,-0.536382e-5};

		if(x==floor(x)&&x<=0) sum=DBL_MAX;					// 0 or negative integer
		else if(x==floor(x))	{											// positive integer
			sum=0;
			for(t=2;t<x-0.1;t+=1.0)	sum+=log(t);	}
		else if(x==0.5)	sum=0.572364942;						// 1/2
		else if(2*x==floor(2*x)&&x>0)	{							// other positive half integer
			sum=0.572364942;
			for(t=0.5;t<x-0.1;t+=1)	sum+=log(t);	}
		else if(2*x==floor(2*x))	{									// negative half integer
			sum=0.572364942;
			for(t=0.5;t<-x+0.1;t+=1)	sum-=log(t);	}
		else if(x<0)																// other negative
			sum=gammalnD(x+1)-log(-x);
		else	{																			// other positive
			x-=1.0;
			t=x+5.5;
			t-=(x+0.5)*log(t);
			sum=1.0;
			for(j=0;j<=5;j++)	{
				x+=1.0;
				sum+=c[j]/x;	}
			sum=-t+log(2.50662827465*sum);	}
		return(sum);	}

	double gammpD(double a,double x)	{ // Returns incomplete gamma function, partly from Numerical Recipes
		double sum,del,ap,eps;
		double gold=0,g=1,fac=1,b1=1,b0=0,anf,ana,an=0,a1,a0=1;

		eps=3e-7;
		if(x<0||a<=0) return -1;			// out of bounds
		else if(x==0) return 0;
		else if(x<a+1)	{
			ap=a;
			del=sum=1.0/a;
			while(fabs(del)>fabs(sum)*eps&&ap-a<100)	{
				ap+=1;
				del*=x/ap;
				sum+=del;	}
			return sum*exp(-x+a*log(x)-gammalnD(a));	}
		else {
			a1=x;
			for(an=1;an<100;an++)	{
				ana=an-a;
				a0=(a1+a0*ana)*fac;
				b0=(b1+b0*ana)*fac;
				anf=an*fac;
				a1=x*a0+anf*a1;
				b1=x*b0+anf*b1;
				if(a1)	{
					fac=1.0/a1;
					g=b1*fac;
					if(fabs((g-gold)/g)<eps)
						return 1.0-exp(-x+a*log(x)-gammalnD(a))*g;
					gold=g;	}}}
		return -1;	}							// shouldn't ever get here

	double erfccD(double x) {
		double t,z,ans;

		z=fabs(x);
		t=1.0/(1.0+0.5*z);
		ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))));
		return x>=0.0?ans:2.0-ans; }

	double erfnD(double x)	{				// Numerical Recipies
		return (x<0?-gammpD(0.5,x*x):gammpD(0.5,x*x));	}

	int iseven(int x)	{
		return(1-abs(x)%2);	}

	void linefitD(double *x,double *y,int n,double *m,double *b) {
		int i;
		double sx,sy,sxx,sxy,denom;

		sx=sy=sxx=sxy=0;
		for(i=0;i<n;i++) {
			sx+=x[i];
			sy+=y[i];
			sxx+=x[i]*x[i];
			sxy+=x[i]*y[i]; }
		denom=n*sxx-sx*sx;
		if(b) *b=(sxx*sy-sx*sxy)/denom;
		if(m) *m=(n*sxy-sx*sy)/denom;
		return; }

double experfcD(double x) {
	double ans,xxinv;
	
	if(fabs(x)<20) ans=exp(x*x)*erfccD(x);
	else {
		xxinv=1.0/(x*x);
		ans=1.0+xxinv*(-1.0/2.0+xxinv*(3.0/4.0+xxinv*(-15.0/8.0+xxinv*(105.0/16.0+xxinv*(-945.0/32.0)))));
		ans/=x*SQRTPI;
		if(x<0) ans+=2.0*exp(x*x); }
	return ans; }

#endif