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

This file calculates reaction rates for a stochastic spatial numerical method,
as described in the research paper "Stochastic Simulation of Chemical Reactions
with Spatial Resolution and Single Molecule Detail" written by Steven S. Andrews
and Dennis Bray, which was published in Physical Biology in 2004.  The code
here is written entirely in ANSII C.

History: Started 3/02, converted to a separate file 5/31/03.
  Several updates made and documentation written 4/08.   */


#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include "rxnparam.h"

#define PI 3.14159265358979323846
#define SQRT2PI 2.50662827462

double rxnparam_erfccD(double x);


/******************************************************************************/
/****************************  UTILITY FUNCTIONS  *****************************/
/******************************************************************************/

/* modelrxnrate */
double modelrxnrate(double a,double b,double difc,double chi) {
	double rate;

	if(b<0) rate=4*PI*difc*a*chi;
	else if(b<=a) rate=-1;
	else rate=4*PI*difc*a*b*chi/(b-a);
	return rate; }


/******************************************************************************/
/***  LOOK-UP FUNCTIONS FOR REACTION RATES AND BINDING AND UNBINDING RADII  ***/
/******************************************************************************/

/* numrxnrate calculates the bimolecular reaction rate for the simulation
algorithm, based on the rms step length in step, the binding radius in a, and
the unbinding radius in b.  It uses cubic polynomial interpolation on previously
computed data, with interpolation both on the reduced step length (step/a) and
on the reduced unbinding radius (b/a).  Enter a negative value of b for
irreversible reactions.  If the input parameters result in reduced values that
are off the edge of the tabulated data, analytical results are used where
possible.  If there are no analytical results, the data are extrapolated.
Variable components: s is step, i is irreversible, r is reversible with b>a, b
is reversible with b<a, y is a rate.  The returned value is the reaction rate
times the time step delta_t.  In other words, this rate constant is per time
step, not per time unit; it has units of length cubed.  */
double numrxnrate(double step,double a,double b) {
	static const double yi[]={				// 600 pts, eps=1e-5, up&down, error<2.4%
		// steps along s in decreasing order
		4.188584,4.188381,4.187971,4.187141,4.185479,4.182243,4.1761835,4.1652925,4.146168,4.112737,4.054173,3.95406,3.7899875,3.536844,3.178271,2.7239775,2.2178055,1.7220855,1.2875175,0.936562,0.668167,0.470021,0.327102,0.225858,0.1551475,0.1061375,0.072355,0.049183,0.0333535,0.022576,0.015257};
	static const double yr[]={				// 500 pts, eps=1e-5, up&down, error<2.2%
		// within a row steps along s in decreasing order; going down rows steps along b in increasing order
		4.188790,4.188790,4.188790,4.188788,4.188785,4.188775,4.188747,4.188665,4.188438,4.187836,4.186351,4.182961,4.175522,4.158666,4.119598,4.036667,3.883389,3.641076,3.316941,2.945750,2.566064,2.203901,1.872874,1.578674,1.322500,1.103003,0.916821,0.759962,0.628537,0.518952,0.428136,
		4.188790,4.188790,4.188789,4.188786,4.188778,4.188756,4.188692,4.188509,4.188004,4.186693,4.183529,4.176437,4.160926,4.125854,4.047256,3.889954,3.621782,3.237397,2.773578,2.289205,1.831668,1.427492,1.087417,0.811891,0.595737,0.430950,0.308026,0.217994,0.152981,0.106569,0.073768,
		4.188790,4.188789,4.188788,4.188783,4.188768,4.188727,4.188608,4.188272,4.187364,4.185053,4.179623,4.167651,4.141470,4.082687,3.956817,3.722364,3.355690,2.877999,2.353690,1.850847,1.411394,1.050895,0.768036,0.553077,0.393496,0.277134,0.193535,0.134203,0.092515,0.063465,0.043358,
		4.188790,4.188789,4.188786,4.188777,4.188753,4.188682,4.188480,4.187919,4.186431,4.182757,4.174363,4.156118,4.116203,4.028391,3.851515,3.546916,3.110183,2.587884,2.055578,1.574898,1.174608,0.858327,0.617223,0.438159,0.307824,0.214440,0.148358,0.102061,0.069885,0.047671,0.032414,
		4.188789,4.188788,4.188783,4.188769,4.188729,4.188614,4.188288,4.187399,4.185108,4.179638,4.167499,4.141379,4.084571,3.964576,3.739474,3.381771,2.906433,2.372992,1.854444,1.401333,1.032632,0.746509,0.531766,0.374442,0.261247,0.180929,0.124557,0.085332,0.058229,0.039604,0.026865,
		4.188789,4.188786,4.188778,4.188756,4.188692,4.188510,4.188002,4.186650,4.183288,4.175566,4.158864,4.123229,4.047257,3.896024,3.632545,3.242239,2.750962,2.220229,1.717239,1.285681,0.939696,0.674540,0.477621,0.334612,0.232460,0.160415,0.110103,0.075242,0.051236,0.034789,0.023565,
		4.188788,4.188784,4.188772,4.188737,4.188637,4.188354,4.187585,4.185607,4.180895,4.170490,4.148460,4.102035,4.006907,3.829897,3.541151,3.133335,2.635739,2.109597,1.619284,1.204298,0.875185,0.625190,0.440882,0.307818,0.213235,0.146800,0.100560,0.068608,0.046655,0.031645,0.021415,
		4.188787,4.188780,4.188761,4.188707,4.188552,4.188124,4.186995,4.184222,4.177931,4.164527,4.136619,4.079198,3.967945,3.772932,3.468711,3.050092,2.548726,2.026932,1.547039,1.144954,0.828592,0.589866,0.414774,0.288895,0.199728,0.137273,0.093903,0.063994,0.043478,0.029467,0.019931,
		4.188785,4.188774,4.188745,4.188661,4.188426,4.187795,4.186203,4.182494,4.174471,4.157859,4.123939,4.056910,3.934028,3.727335,3.412145,2.985387,2.481620,1.963935,1.492519,1.100533,0.793976,0.563797,0.395615,0.275071,0.189896,0.130359,0.089086,0.060661,0.041187,0.027899,0.018863,
		4.188782,4.188766,4.188720,4.188592,4.188244,4.187347,4.185205,4.180486,4.170691,4.150875,4.111574,4.037574,3.906673,3.691217,3.367270,2.934386,2.429282,1.915212,1.450660,1.066615,0.767733,0.544143,0.381225,0.264724,0.182558,0.125209,0.085504,0.058187,0.039490,0.026741,0.018074,
		4.188777,4.188752,4.188683,4.188492,4.187993,4.186777,4.184049,4.178342,4.166861,4.144167,4.100756,4.021817,3.884852,3.662169,3.331386,2.893934,2.388069,1.877097,1.418051,1.040351,0.747549,0.529083,0.370239,0.256844,0.176983,0.121304,0.082791,0.056318,0.038211,0.025871,0.017486,
		4.188770,4.188732,4.188628,4.188352,4.187671,4.186115,4.182830,4.176234,4.163261,4.138284,4.092055,4.009206,3.867171,3.638731,3.302580,2.861651,2.355385,1.846979,1.392364,1.019841,0.731849,0.517409,0.361743,0.250766,0.172684,0.118295,0.080706,0.054883,0.037230,0.025204,0.017035,
		4.188759,4.188702,4.188551,4.188171,4.187294,4.185421,4.181657,4.174292,4.160089,4.133434,4.084980,3.998932,3.852801,3.619743,3.279339,2.835771,2.329266,1.822929,1.372041,1.003734,0.719559,0.508296,0.355130,0.246037,0.169346,0.115966,0.079096,0.053783,0.036482,0.024699,0.016694,
		4.188742,4.188659,4.188449,4.187958,4.186900,4.184765,4.180606,4.172596,4.157440,4.129551,4.079183,3.990499,3.841031,3.604270,3.260519,2.814871,2.308170,1.803653,1.355971,0.991032,0.709899,0.501150,0.349948,0.242337,0.166741,0.114156,0.077855,0.052940,0.035912,0.024316,0.016437,
		4.188719,4.188603,4.188330,4.187736,4.186534,4.184200,4.179711,4.171172,4.155254,4.126372,4.074457,3.983648,3.831433,3.591609,3.245083,2.797662,2.290945,1.788288,1.343237,0.981003,0.702295,0.495533,0.345879,0.239437,0.164709,0.112756,0.076902,0.052295,0.035478,0.024023,0.016239,
		4.188688,4.188536,4.188204,4.187532,4.186227,4.183725,4.178948,4.169962,4.153461,4.123737,4.070508,3.977891,3.823388,3.580958,3.232035,2.783281,2.277011,1.776032,1.333144,0.973094,0.696303,0.491107,0.342676,0.237173,0.163143,0.111689,0.076180,0.051808,0.035151,0.023802,0.016089};
	const double yb[]={							// 100 pts, eps=1e-5, down only
		// within a row steps along s in decreasing order; going down rows steps along b in increasing order
		4.188790,4.188791,4.188791,4.188793,4.188798,4.188814,4.188859,4.188986,4.189328,4.190189,4.192213,4.196874,4.208210,4.236983,4.307724,4.473512,4.852754,5.723574,7.838600,13.880005,40.362527,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		4.188790,4.188791,4.188791,4.188793,4.188798,4.188813,4.188858,4.188982,4.189319,4.190165,4.192156,4.196738,4.207879,4.236133,4.305531,4.467899,4.838139,5.683061,7.710538,13.356740,36.482501,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		4.188790,4.188791,4.188791,4.188793,4.188798,4.188812,4.188854,4.188973,4.189292,4.190095,4.191984,4.196332,4.206887,4.233594,4.298997,4.451212,4.795122,5.566006,7.352764,11.994360,28.085852,261.834153,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		4.188790,4.188791,4.188791,4.188792,4.188797,4.188810,4.188848,4.188956,4.189246,4.189978,4.191699,4.195657,4.205243,4.229397,4.288214,4.424093,4.726265,5.384525,6.831926,10.241985,19.932261,69.581522,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		4.188790,4.188791,4.188791,4.188792,4.188796,4.188807,4.188840,4.188933,4.189183,4.189814,4.191300,4.194716,4.202957,4.223593,4.273478,4.387310,4.635278,5.155725,6.228070,8.494834,13.840039,30.707995,217.685066,-1,-1,-1,-1,-1,-1,-1,-1,
		4.188790,4.188791,4.188791,4.188792,4.188795,4.188803,4.188829,4.188903,4.189101,4.189604,4.190789,4.193514,4.200048,4.216255,4.254956,4.342131,4.526936,4.897901,5.610205,6.966144,9.704274,16.182908,37.778925,387.192994,-1,-1,-1,-1,-1,-1,-1,
		4.188790,4.188790,4.188791,4.188791,4.188793,4.188799,4.188817,4.188866,4.189002,4.189348,4.190168,4.192054,4.196535,4.207468,4.233117,4.289818,4.406055,4.627980,5.025468,5.719284,6.973164,9.454673,15.165237,32.669841,175.505441,-1,-1,-1,-1,-1,-1,
		4.188790,4.188790,4.188790,4.188791,4.188791,4.188794,4.188801,4.188823,4.188885,4.189047,4.189439,4.190344,4.192444,4.197336,4.208277,4.231803,4.277545,4.359611,4.499274,4.738728,5.168746,5.973516,7.554064,10.974866,20.015747,59.729701,-1,-1,-1,-1,-1,
		4.188790,4.188790,4.188790,4.188790,4.188789,4.188788,4.188784,4.188774,4.188751,4.188701,4.188602,4.188390,4.187803,4.185972,4.180908,4.169592,4.145768,4.102670,4.041148,3.980771,3.961901,4.039620,4.292076,4.858676,6.044956,8.672140,15.704344,48.437802,-1,-1,-1,
		4.188790,4.188790,4.188790,4.188789,4.188787,4.188781,4.188764,4.188718,4.188599,4.188311,4.187661,4.186201,4.182644,4.173501,4.151495,4.104610,4.014334,3.863341,3.650600,3.398722,3.141206,2.906174,2.713737,2.580237,2.524103,2.574058,2.785372,3.278275,4.350631,6.897568,14.924070,
		4.188790,4.188790,4.188790,4.188788,4.188784,4.188773,4.188742,4.188656,4.188430,4.187878,4.186617,4.183786,4.177000,4.160053,4.120473,4.038061,3.886243,3.645077,3.322103,2.951900,2.573084,2.212134,1.883264,1.592189,1.339766,1.124235,0.942582,0.791354,0.667141,0.566827,0.487862};
	const double slo=3,sinc=-0.2;								// logs of values; shi=-3
	const double blor=0,bincr=0.2;								// logs of values; bhir=3
	const double blob=0,bincb=0.1;								// actual values; bhib=1
	const int snum=31,bnumr=16,bnumb=11;
	double x[4],y[4],z[4];
	int sindx,bindx,i,j;
	double ans;

	if(step<0 || a<0) return -1;
	if(step==0 && b>=0 && b<1) return -1;
	if(step==0) return 0;
	step=log(step/a);
	b/=a;

	sindx=(int)floor((step-slo)/sinc);
	for(i=0;i<4;i++) x[i]=slo+(sindx-1+i)*sinc;
	z[0]=(step-x[1])*(step-x[2])*(step-x[3])/(-6.0*sinc*sinc*sinc);
	z[1]=(step-x[0])*(step-x[2])*(step-x[3])/(2.0*sinc*sinc*sinc);
	z[2]=(step-x[0])*(step-x[1])*(step-x[3])/(-2.0*sinc*sinc*sinc);
	z[3]=(step-x[0])*(step-x[1])*(step-x[2])/(6.0*sinc*sinc*sinc);

	if(b<0)
		for(ans=i=0;i<4;i++) {
			if(sindx-1+i>=snum) ans+=z[i]*2.0*PI*exp(2.0*x[i]);
			else if(sindx-1+i<0) ans+=z[i]*4.0*PI/3.0;
			else ans+=z[i]*yi[sindx-1+i]; }
	else if(b<1.0) {
		bindx=(int)floor((b-blob)/bincb);
		if(bindx<1) bindx=1;
		else if(bindx+2>=bnumb) bindx=bnumb-3;
		while(sindx+3>=snum || (sindx>0 && yb[(bindx-1)*snum+sindx+3]<0)) sindx--;
		for(i=0;i<4;i++) x[i]=slo+(sindx-1+i)*sinc;
		z[0]=(step-x[1])*(step-x[2])*(step-x[3])/(-6.0*sinc*sinc*sinc);
		z[1]=(step-x[0])*(step-x[2])*(step-x[3])/(2.0*sinc*sinc*sinc);
		z[2]=(step-x[0])*(step-x[1])*(step-x[3])/(-2.0*sinc*sinc*sinc);
		z[3]=(step-x[0])*(step-x[1])*(step-x[2])/(6.0*sinc*sinc*sinc);
		for(j=0;j<4;j++)
			for(y[j]=i=0;i<4;i++) {
				if(sindx-1+i>=snum) y[j]+=z[i]*yb[(bindx-1+j)*snum];
				else if(sindx-1+i<0) y[j]+=z[i]*4.0*PI/3.0;
				else y[j]+=z[i]*yb[(bindx-1+j)*snum+sindx-1+i]; }
		for(j=0;j<4;j++) x[j]=blob+(bindx-1+j)*bincb;
		z[0]=(b-x[1])*(b-x[2])*(b-x[3])/(-6.0*bincb*bincb*bincb);
		z[1]=(b-x[0])*(b-x[2])*(b-x[3])/(2.0*bincb*bincb*bincb);
		z[2]=(b-x[0])*(b-x[1])*(b-x[3])/(-2.0*bincb*bincb*bincb);
		z[3]=(b-x[0])*(b-x[1])*(b-x[2])/(6.0*bincb*bincb*bincb);
		ans=z[0]*y[0]+z[1]*y[1]+z[2]*y[2]+z[3]*y[3]; }
	else {
		b=log(b);
		bindx=(int)floor((b-blor)/bincr);
		if(bindx<1) bindx=1;
		else if(bindx+2>=bnumr) bindx=bnumr-3;
		for(j=0;j<4;j++)
			for(y[j]=i=0;i<4;i++) {
				if(sindx-1+i>=snum && b==0) y[j]+=z[i]*2.0*PI*exp(x[i])*(1.0+exp(x[i]));
				else if(sindx-1+i>=snum) y[j]+=z[i]*2.0*PI*exp(2.0*x[i])*exp(b)/(exp(b)-1.0);
				else if(sindx-1+i<0) y[j]+=z[i]*4.0*PI/3.0;
				else y[j]+=z[i]*yr[(bindx-1+j)*snum+sindx-1+i]; }
		for(j=0;j<4;j++) x[j]=blor+(bindx-1+j)*bincr;
		z[0]=(b-x[1])*(b-x[2])*(b-x[3])/(-6.0*bincr*bincr*bincr);
		z[1]=(b-x[0])*(b-x[2])*(b-x[3])/(2.0*bincr*bincr*bincr);
		z[2]=(b-x[0])*(b-x[1])*(b-x[3])/(-2.0*bincr*bincr*bincr);
		z[3]=(b-x[0])*(b-x[1])*(b-x[2])/(6.0*bincr*bincr*bincr);
		ans=z[0]*y[0]+z[1]*y[1]+z[2]*y[2]+z[3]*y[3]; }
	return ans*a*a*a; }



/* numrxnrateprob */
double numrxnrateprob(double step,double a,double b,double prob) {
	static const double kirr[]={				// 200 pts, eps=1e-6, up&down, error<4.9%, av.error 0.36%
		// within a row steps along s in increasing order; going down rows steps along p in increasing order
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.010134,0.0134585,0.0172905,0.0214105,0.0255145,0.029298,0.0325455,0.035164,0.0371705,0.038649,0.039706,0.0404445,0.0409505,0.041292,0.041518,0.0416635,0.041755,0.0418105,0.041844,0.041863,0.041874,0.04188,0.041884,0.041886,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0117075,0.016229,0.0219785,0.0289295,0.036822,0.0451515,0.0532835,0.060637,0.0668365,0.07176,0.075484,0.0781955,0.0801095,0.081426,0.082309,0.082884,0.0832465,0.083468,0.0836,0.0836765,0.08372,0.0837445,0.083759,0.083766,0.083771,0.083773,0.083774,0.083775,0.083775,0.083776,0.083776,
		0.013152,0.0188265,0.0266425,0.0371035,0.0506205,0.0672785,0.0866085,0.1075055,0.128414,0.1477555,0.1643765,0.177756,0.1879465,0.195351,0.200508,0.2039535,0.2061635,0.2075285,0.2083445,0.2088205,0.2090925,0.2092465,0.209332,0.2093805,0.209408,0.209423,0.209431,0.209435,0.209437,0.209438,0.209439,
		0.013906,0.020166,0.0290685,0.0415065,0.058521,0.081166,0.1101695,0.145482,0.185849,0.2287295,0.2708075,0.3089305,0.340963,0.3660955,0.384619,0.397497,0.405978,0.411304,0.414522,0.4164105,0.417495,0.418107,0.4184495,0.418643,0.418753,0.418814,0.418846,0.418863,0.418871,0.418875,0.418877,
		0.014458,0.0211455,0.0308165,0.0446825,0.0642925,0.0916025,0.12884,0.1781235,0.2407515,0.316199,0.4012615,0.490054,0.5752905,0.6503725,0.7111165,0.756394,0.7876975,0.8079965,0.820509,0.8279405,0.8322395,0.834675,0.836042,0.836815,0.837253,0.837496,0.837626,0.837692,0.837725,0.837742,0.83775,
		0.0149755,0.022065,0.032443,0.0475935,0.0695775,0.101199,0.146253,0.2096285,0.2971995,0.415231,0.5688585,0.75941,0.9810695,1.218721,1.449651,1.6506675,1.807287,1.9178015,1.98982,2.0340965,2.0602395,2.0752335,2.083706,2.088511,2.0912435,2.09276,2.093569,2.093984,2.094191,2.094293,2.094344,
		0.015258,0.0225675,0.0333325,0.049161,0.0723835,0.1062885,0.155483,0.226384,0.3276425,0.4704975,0.6686825,0.937162,1.2882245,1.722884,2.2186495,2.724778,3.178969,3.5373775,3.790392,3.9543465,4.054361,4.112839,4.1462465,4.1653235,4.1762035,4.182256,4.185487,4.187145,4.187973,4.188382,4.188584};
	static const double krev[]={				// 200 pts, eps=1e-6, up&down, error<5.3%, averror 0.21%
		// within a row steps along s in increasing order; going down rows steps along p in increasing order; going down groups steps along b in increasing order
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.028984,0.0318765,0.0343805,0.0364265,0.038016,0.039201,0.040057,0.040661,0.0410795,0.041365,0.041558,0.041686,0.0417695,0.041823,0.0418545,0.041872,0.041881,0.041885,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0470995,0.0535025,0.059681,0.065288,0.07007,0.073924,0.076882,0.0790655,0.080629,0.081722,0.0824715,0.082976,0.083307,0.083517,0.083643,0.0837125,0.083748,0.083764,0.08377,0.083773,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.083615,0.097945,0.1134355,0.1294995,0.145305,0.159944,0.1726755,0.1831025,0.1912045,0.1972265,0.2015405,0.204533,0.2065415,0.2078335,0.208613,0.209046,0.2092645,0.209364,0.209406,0.209425,0.209433,0.209437,0.209439,0.209439,0.209439,0.209439,0.209439,0.20944,0.20944,0.20944,0.20944,
		0.1250115,0.14843,0.175039,0.20455,0.2362085,0.2687045,0.300276,0.3290955,0.3537695,0.37365,0.388818,0.3998415,0.407488,0.4125145,0.4155885,0.4173085,0.4181805,0.418577,0.418746,0.418821,0.418855,0.418869,0.418875,0.418878,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,
		0.1839805,0.2202195,0.262607,0.311499,0.3668445,0.4279295,0.49304,0.559275,0.6228695,0.680067,0.7281325,0.7659105,0.793713,0.8127615,0.824716,0.831505,0.8349685,0.8365495,0.837226,0.837524,0.83766,0.83772,0.837744,0.837753,0.837756,0.837757,0.837758,0.837758,0.837758,0.837758,0.837758,
		0.3023935,0.3639145,0.4375255,0.52495,0.627794,0.7473325,0.884173,1.037587,1.2046405,1.379149,1.5515315,1.7101205,1.8440585,1.946143,2.0149835,2.055779,2.0770625,2.086855,2.0910715,2.092935,2.093784,2.094156,2.094307,2.094364,2.094384,2.094391,2.094394,2.094395,2.094395,2.094395,2.094395,
		0.438632,0.5281975,0.6367805,0.767722,0.924435,1.1105185,1.3296165,1.584829,1.877914,2.2079835,2.5693915,2.9484465,3.319073,3.6427,3.884565,4.037486,4.120112,4.1587575,4.175513,4.182953,4.186348,4.187835,4.188438,4.188665,4.188746,4.188775,4.188785,4.188788,4.18879,4.18879,4.18879,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.021711,0.0255795,0.0292035,0.032373,0.0349775,0.037006,0.0385195,0.039613,0.040381,0.040911,0.041269,0.0415075,0.041663,0.041762,0.0418215,0.0418555,0.041873,0.041882,0.041885,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.030514,0.0378765,0.045654,0.053337,0.0604115,0.0664965,0.0714195,0.0752005,0.0779855,0.0799685,0.081342,0.082271,0.082883,0.083274,0.083512,0.083647,0.083717,0.08375,0.083765,0.083771,0.083774,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.042629,0.0558905,0.071724,0.089781,0.109258,0.128942,0.147471,0.1637125,0.177035,0.187346,0.1949425,0.2003095,0.203957,0.206335,0.2078,0.2086375,0.209074,0.20928,0.209369,0.209408,0.209426,0.209434,0.209438,0.209439,0.209439,0.209439,0.209439,0.20944,0.20944,0.20944,0.20944,
		0.0514115,0.069516,0.0924725,0.120699,0.1541035,0.1917935,0.231883,0.2716955,0.308419,0.3398965,0.365095,0.384065,0.3975705,0.4066585,0.4123755,0.4156845,0.417422,0.4182415,0.418598,0.418755,0.418826,0.418858,0.418871,0.418876,0.418878,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,
		0.059424,0.082353,0.1126975,0.1520045,0.201619,0.26228,0.333503,0.4129205,0.4960515,0.576961,0.649794,0.710392,0.757085,0.7903865,0.8121685,0.825082,0.8319545,0.8352115,0.836635,0.837261,0.837547,0.837674,0.837727,0.837747,0.837754,0.837757,0.837758,0.837758,0.837758,0.837758,0.837758,
		0.068507,0.097392,0.1372535,0.1914965,0.264124,0.3594505,0.481498,0.633012,0.8140055,1.0200345,1.240778,1.460258,1.6598735,1.823422,1.941938,2.0170155,2.0585055,2.0785525,2.0873905,2.0912925,2.093076,2.09387,2.094199,2.094325,2.09437,2.094387,2.094392,2.094394,2.094395,2.094395,2.094395,
		0.074197,0.1070925,0.153612,0.218761,0.3089765,0.4321595,0.5973065,0.8137355,1.089582,1.4299205,1.8340335,2.291319,2.7753935,3.2388975,3.6229535,3.890774,4.0477635,4.1259225,4.1608605,4.176397,4.1835175,4.186691,4.188005,4.188509,4.188692,4.188756,4.188778,4.188786,4.188789,4.18879,4.18879,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.017987,0.021987,0.0259585,0.029623,0.032773,0.035321,0.0372785,0.038725,0.039762,0.040488,0.040988,0.0413255,0.0415495,0.0416935,0.041783,0.0418355,0.041863,0.041877,0.041883,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.023636,0.030497,0.0381905,0.046258,0.0541225,0.0612435,0.0672655,0.072064,0.0757075,0.07837,0.0802575,0.0815595,0.0824335,0.0830035,0.083358,0.0835655,0.0836775,0.083732,0.083757,0.083767,0.083772,0.083774,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.030314,0.041186,0.054875,0.0713915,0.090284,0.110555,0.130793,0.1495435,0.1657185,0.178807,0.1888295,0.196146,0.201261,0.2046815,0.2068505,0.2081325,0.2088275,0.2091675,0.20932,0.209386,0.209416,0.20943,0.209436,0.209438,0.209439,0.209439,0.209439,0.209439,0.20944,0.20944,0.20944,
		0.03451,0.0481425,0.066248,0.089653,0.118868,0.153756,0.193154,0.234782,0.275643,0.312821,0.3442505,0.3690655,0.3874615,0.40029,0.4086535,0.413685,0.41644,0.417794,0.4184,0.418666,0.418787,0.418842,0.418865,0.418874,0.418877,0.418878,0.418879,0.418879,0.418879,0.418879,0.418879,
		0.037948,0.0539735,0.0760255,0.1058445,0.145278,0.195976,0.258855,0.333305,0.4164795,0.5031695,0.5867225,0.660845,0.7213265,0.766712,0.7978675,0.8172425,0.828062,0.833428,0.835845,0.8369065,0.83739,0.837608,0.837701,0.837737,0.837751,0.837756,0.837757,0.837758,0.837758,0.837758,0.837758,
		0.041462,0.060056,0.0864675,0.123597,0.17514,0.2455735,0.3399605,0.463436,0.6200605,0.8107845,1.0306265,1.2665115,1.498095,1.70265,1.8622995,1.9708425,2.034847,2.0675545,2.082479,2.089082,2.0920955,2.093459,2.094038,2.094266,2.09435,2.094379,2.09439,2.094393,2.094394,2.094395,2.094395,
		0.0434845,0.063618,0.0927,0.1344305,0.193821,0.2775065,0.393983,0.5536985,0.768791,1.0517725,1.4123445,1.8518005,2.354633,2.878803,3.356234,3.722846,3.9572635,4.082812,4.141392,4.1675915,4.1796015,4.185049,4.187364,4.188273,4.188608,4.188727,4.188768,4.188783,4.188788,4.188789,4.18879,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0157705,0.0197185,0.023794,0.0276955,0.0311645,0.034051,0.0363205,0.0380265,0.0392655,0.040141,0.0407495,0.041164,0.0414425,0.041626,0.041742,0.0418115,0.041851,0.041871,0.041881,0.041885,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0199525,0.026301,0.033682,0.041724,0.049872,0.0575245,0.064209,0.069682,0.073927,0.0770805,0.079346,0.0809315,0.082014,0.082734,0.0831945,0.083473,0.083629,0.083709,0.083746,0.083763,0.08377,0.083773,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.024516,0.033887,0.0460225,0.061138,0.0790465,0.0990015,0.119714,0.1396385,0.157419,0.172233,0.1838635,0.192555,0.1987775,0.203051,0.205844,0.207557,0.2085265,0.2090235,0.209255,0.209357,0.209403,0.209424,0.209434,0.209437,0.209439,0.209439,0.209439,0.209439,0.20944,0.20944,0.20944,
		0.027193,0.0384625,0.053765,0.0740565,0.1001285,0.132287,0.169931,0.2112565,0.253421,0.2932375,0.3280955,0.3565555,0.3783665,0.394103,0.404745,0.4114175,0.4152435,0.4172185,0.4181405,0.418549,0.418734,0.418819,0.418855,0.41887,0.418876,0.418878,0.418879,0.418879,0.418879,0.418879,0.418879,
		0.029287,0.042099,0.060034,0.08477,0.1182345,0.1623885,0.218787,0.287807,0.3677535,0.4543585,0.541299,0.621796,0.6904505,0.744344,0.783111,0.808443,0.823342,0.831146,0.834809,0.8364385,0.8371775,0.837516,0.837664,0.837723,0.837746,0.837754,0.837757,0.837758,0.837758,0.837758,0.837758,
		0.0313395,0.0457155,0.066369,0.0957945,0.137288,0.1950325,0.274055,0.3799305,0.517891,0.691126,0.898245,1.130586,1.37101,1.596335,1.783943,1.920464,2.00659,2.0535645,2.076061,2.086168,2.0907705,2.092883,2.093805,2.094177,2.094318,2.094368,2.094386,2.094392,2.094394,2.094395,2.094395,
		0.032484,0.0477545,0.0699855,0.1021835,0.1485155,0.2146545,0.308122,0.4385455,0.617689,0.8588385,1.175085,1.5753325,2.056003,2.5882955,3.11056,3.547265,3.8518625,4.0285985,4.116088,4.1560095,4.174317,4.182746,4.186431,4.18792,4.188481,4.188682,4.188753,4.188777,4.188786,4.188789,4.18879,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0143235,0.0181815,0.0222725,0.0262945,0.029961,0.0330775,0.0355725,0.037473,0.038868,0.039862,0.0405555,0.0410315,0.041352,0.041565,0.0417025,0.041788,0.0418375,0.041864,0.041877,0.041883,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0176945,0.0236365,0.0307125,0.0386245,0.0468595,0.054801,0.0619065,0.0678465,0.072531,0.0760575,0.0786165,0.0804195,0.081661,0.082495,0.0830385,0.083377,0.083575,0.083681,0.083733,0.083757,0.083767,0.083772,0.083774,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.021198,0.029593,0.040653,0.0547065,0.071737,0.0912,0.111951,0.1324565,0.1512235,0.1672085,0.1799915,0.1896845,0.196715,0.201617,0.204892,0.206965,0.2081885,0.2088485,0.2091725,0.209319,0.209386,0.209417,0.20943,0.209436,0.209438,0.209439,0.209439,0.209439,0.209439,0.20944,0.20944,
		0.0231735,0.0330275,0.0465805,0.0648255,0.088684,0.1187185,0.1547045,0.195242,0.23774,0.278966,0.3159655,0.346836,0.370962,0.3887365,0.401081,0.409098,0.413905,0.416525,0.4178115,0.4183985,0.418664,0.418787,0.418842,0.418865,0.418874,0.418877,0.418878,0.418879,0.418879,0.418879,0.418879,
		0.0246795,0.0356775,0.0512175,0.0728905,0.102602,0.142411,0.1941835,0.2588805,0.33563,0.4209885,0.509057,0.592823,0.6661845,0.725433,0.7695105,0.7995365,0.818098,0.828396,0.833499,0.8358375,0.8368995,0.83739,0.83761,0.837702,0.837738,0.837751,0.837756,0.837757,0.837758,0.837758,0.837758,
		0.0261265,0.0382475,0.055766,0.0809045,0.116656,0.1669175,0.236524,0.3311005,0.456386,0.616768,0.8128145,1.038307,1.2785415,1.51185,1.7149385,1.8709735,1.9757305,2.0368515,2.0679825,2.0824335,2.089037,2.092098,2.093473,2.094047,2.09427,2.094351,2.09438,2.09439,2.094393,2.094394,2.094395,
		0.0269215,0.03967,0.058305,0.0854235,0.124673,0.181092,0.261484,0.3747615,0.532164,0.746946,1.0330375,1.401716,1.8548455,2.373404,2.9068175,3.3821345,3.7397705,3.9648045,4.084467,4.1412135,4.1674115,4.1796125,4.185104,4.1874,4.188289,4.188614,4.188729,4.188769,4.188783,4.188788,4.188789,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0133205,0.017089,0.021163,0.0252485,0.0290425,0.032321,0.0349825,0.0370325,0.038549,0.0396365,0.040399,0.0409235,0.0412785,0.0415155,0.041669,0.041766,0.0418235,0.0418565,0.041873,0.041881,0.041885,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.016194,0.021825,0.0286435,0.03641,0.044651,0.0527555,0.060141,0.066414,0.0714275,0.07524,0.078029,0.080007,0.081375,0.082298,0.0829045,0.0832885,0.0835195,0.0836495,0.083717,0.083749,0.083764,0.08377,0.083774,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.019085,0.0268125,0.0371075,0.050367,0.0666895,0.085674,0.106308,0.1271055,0.146503,0.1633085,0.1769425,0.187403,0.195063,0.200446,0.204079,0.2064215,0.207848,0.2086555,0.209074,0.209273,0.209364,0.209406,0.209426,0.209434,0.209438,0.209439,0.209439,0.209439,0.209439,0.20944,0.20944,
		0.020674,0.0296055,0.041988,0.058822,0.081097,0.109524,0.1441335,0.183835,0.226279,0.268277,0.3066895,0.339285,0.3651295,0.384406,0.397978,0.4069785,0.412561,0.4157565,0.4174195,0.4182135,0.418577,0.418746,0.418824,0.418858,0.418871,0.418876,0.418878,0.418879,0.418879,0.418879,0.418879,
		0.021868,0.0317215,0.045725,0.0653915,0.092583,0.129385,0.177817,0.239203,0.3132345,0.397114,0.4854025,0.5710995,0.647607,0.7104965,0.758169,0.791482,0.8128625,0.82536,0.831942,0.8351005,0.8365525,0.837227,0.837538,0.837672,0.837727,0.837747,0.837754,0.837757,0.837758,0.837758,0.837758,
		0.0230035,0.0337455,0.0493255,0.071781,0.103884,0.149306,0.2126905,0.2995905,0.415955,0.566848,0.754141,0.9734555,1.211822,1.448389,1.659613,1.827461,1.945473,2.0185935,2.0584295,2.0778625,2.0868755,2.091079,2.093018,2.09386,2.094198,2.094325,2.094371,2.094387,2.094392,2.094394,2.094395,
		0.023626,0.034857,0.051312,0.075328,0.1102085,0.1605585,0.232672,0.3349065,0.478002,0.6749715,0.940094,1.286062,1.717642,2.2206485,2.7513725,3.2426025,3.6328475,3.896208,4.0473645,4.123173,4.158818,4.1755465,4.183284,4.186649,4.188002,4.18851,4.188692,4.188756,4.188778,4.188786,4.188789,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0126,0.016287,0.020334,0.024452,0.0283315,0.031728,0.0345135,0.036679,0.038291,0.0394535,0.0402715,0.0408355,0.0412185,0.0414745,0.041641,0.041747,0.041811,0.0418485,0.041869,0.041879,0.041884,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.015146,0.020538,0.027147,0.034778,0.042994,0.051193,0.05877,0.065286,0.070549,0.0745845,0.077555,0.0796725,0.0811425,0.082139,0.082796,0.0832145,0.0834705,0.083619,0.0837,0.08374,0.083759,0.083768,0.083773,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.017649,0.0249015,0.0346375,0.0472985,0.0630605,0.081628,0.1020985,0.1230385,0.142854,0.16025,0.174523,0.1855765,0.1937325,0.1995015,0.20342,0.205966,0.2075415,0.208463,0.208965,0.209218,0.209338,0.209393,0.20942,0.209432,0.209436,0.209438,0.209439,0.209439,0.209439,0.209439,0.20944,
		0.0190005,0.027295,0.038857,0.05468,0.075793,0.102999,0.136503,0.1754485,0.2176915,0.260122,0.299493,0.3333455,0.360496,0.38095,0.3954785,0.405212,0.4113555,0.4149915,0.416987,0.417994,0.4184715,0.418695,0.4188,0.418847,0.418867,0.418875,0.418877,0.418878,0.418879,0.418879,0.418879,
		0.020003,0.029084,0.0420385,0.060316,0.085733,0.1203775,0.166348,0.225198,0.2970175,0.379502,0.46762,0.554469,0.633173,0.698779,0.7491505,0.7848275,0.808197,0.8223525,0.830228,0.8342265,0.836129,0.837022,0.837441,0.837631,0.83771,0.837741,0.837752,0.837756,0.837757,0.837758,0.837758,
		0.020944,0.0307725,0.0450615,0.065713,0.0953395,0.1374375,0.1964885,0.277946,0.3878395,0.531635,0.7120735,0.9261115,1.1622465,1.4005145,1.6170015,1.7923745,1.918962,2.000701,2.0479555,2.072467,2.0842455,2.089801,2.092417,2.093598,2.094094,2.094286,2.094357,2.094382,2.09439,2.094394,2.094395,
		0.021451,0.031688,0.046708,0.068671,0.10064,0.146915,0.213419,0.308089,0.4412485,0.625636,0.8756175,1.204704,1.619721,2.1100605,2.636198,3.133754,3.54149,3.8301035,4.007068,4.101986,4.1483855,4.170453,4.1808845,4.185605,4.187585,4.188354,4.188637,4.188737,4.188772,4.188784,4.188788,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.012066,0.015683,0.019699,0.0238345,0.0277745,0.0312575,0.034139,0.0363945,0.0380835,0.0393055,0.0401675,0.040764,0.04117,0.0414415,0.0416185,0.041732,0.041801,0.0418415,0.0418645,0.041877,0.041883,0.041885,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0143855,0.019593,0.026031,0.0335455,0.0417255,0.049981,0.057694,0.0643915,0.069846,0.074056,0.077171,0.0794005,0.0809535,0.0820085,0.082707,0.0831545,0.0834295,0.0835915,0.0836825,0.083731,0.083755,0.083766,0.083771,0.083774,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.016628,0.023531,0.032849,0.04505,0.0603705,0.078591,0.0988935,0.1198995,0.140001,0.1578305,0.1725915,0.184108,0.192657,0.1987365,0.202886,0.2055975,0.207288,0.208292,0.2088575,0.2091585,0.209307,0.209378,0.209412,0.209428,0.209435,0.209438,0.209439,0.209439,0.209439,0.209439,0.20944,
		0.0178235,0.02566,0.036625,0.051702,0.071941,0.098212,0.1308365,0.169135,0.2111355,0.253808,0.2938505,0.3286375,0.3567905,0.3781685,0.3934655,0.403789,0.4103615,0.4143145,0.416561,0.4177555,0.4183495,0.418634,0.41877,0.418833,0.418861,0.418872,0.418877,0.418878,0.418879,0.418879,0.418879,
		0.0187045,0.0272375,0.039443,0.05672,0.0808415,0.1138915,0.158011,0.2149045,0.2849475,0.366213,0.45401,0.5415645,0.621831,0.6894765,0.7419585,0.779507,0.804368,0.8196965,0.8285365,0.833278,0.835644,0.8367795,0.837321,0.837575,0.837686,0.837731,0.837749,0.837755,0.837757,0.837758,0.837758,
		0.0195275,0.0287175,0.0420985,0.0614765,0.0893415,0.1290535,0.1849705,0.262442,0.3675215,0.5059255,0.680995,0.8906685,1.1245965,1.363639,1.5838645,1.7648605,1.8975105,1.985054,2.037698,2.0666215,2.081235,2.0882885,2.0916685,2.093251,2.093946,2.094229,2.094336,2.094374,2.094388,2.094393,2.094394,
		0.0199685,0.0295155,0.0435365,0.064064,0.0939895,0.137388,0.1999055,0.2891665,0.415152,0.590347,0.8291035,1.145424,1.5475265,2.0274625,2.5492485,3.0506085,3.4691175,3.773164,3.9679865,4.0791495,4.1364745,4.1644355,4.1778965,4.184215,4.186995,4.188124,4.188552,4.188707,4.188761,4.18878,4.188787,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0116635,0.015221,0.0192055,0.023349,0.027332,0.0308815,0.0338385,0.036165,0.037914,0.0391845,0.040083,0.0407055,0.04113,0.0414145,0.0416005,0.0417195,0.0417925,0.041836,0.0418605,0.041874,0.041881,0.041885,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.01382,0.0188825,0.0251835,0.0325965,0.04074,0.0490305,0.056842,0.0636775,0.069281,0.0736295,0.07686,0.079179,0.080799,0.0819025,0.0826345,0.083105,0.083396,0.083569,0.083667,0.0837215,0.083749,0.083763,0.08377,0.083773,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.015879,0.02252,0.03152,0.043363,0.0583325,0.07627,0.0964185,0.1174485,0.1377505,0.155906,0.1710435,0.182924,0.1917865,0.198114,0.202451,0.205298,0.2070835,0.208152,0.2087635,0.2090995,0.2092745,0.209361,0.209403,0.209424,0.209433,0.209437,0.209439,0.209439,0.209439,0.209439,0.209439,
		0.0169675,0.0244655,0.0349865,0.0495005,0.069069,0.094615,0.126541,0.1643,0.2060595,0.248867,0.2893915,0.324884,0.353816,0.375923,0.3918335,0.4026355,0.4095595,0.413759,0.416185,0.4175215,0.41822,0.4185665,0.418735,0.418816,0.418853,0.418869,0.418875,0.418878,0.418879,0.418879,0.418879,
		0.017768,0.0259,0.037554,0.05409,0.077242,0.1090855,0.1517935,0.207164,0.2757835,0.3560175,0.4434545,0.5314485,0.612852,0.6820495,0.7361765,0.77522,0.801293,0.817527,0.827049,0.8323475,0.8351265,0.836509,0.8371815,0.837505,0.837654,0.837718,0.837743,0.837753,0.837756,0.837757,0.837758,
		0.018515,0.027243,0.0399645,0.0584115,0.0849845,0.122933,0.1765165,0.251004,0.352432,0.4866845,0.657528,0.8636385,1.0955705,1.334891,1.557747,1.7430385,1.8804905,1.972375,2.028727,2.060905,2.0780265,2.086605,2.090795,2.092815,2.093747,2.094146,2.094304,2.094363,2.094384,2.094391,2.094394,
		0.0189185,0.027969,0.041269,0.0607565,0.089196,0.1304945,0.1900855,0.275357,0.3960195,0.5643335,0.794597,1.10112,1.4931035,1.964564,2.4822685,2.985991,3.412643,3.7276185,3.934113,4.056943,4.1238235,4.1577435,4.1744135,4.182476,4.1862,4.187795,4.188426,4.188661,4.188745,4.188774,4.188785,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0113565,0.0148645,0.0188195,0.0229625,0.0269765,0.0305785,0.033594,0.035978,0.0377765,0.039086,0.040014,0.040658,0.0410975,0.041392,0.0415855,0.0417095,0.0417855,0.041831,0.0418575,0.041872,0.04188,0.041884,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.013393,0.018343,0.0245325,0.031857,0.0399645,0.0482785,0.056163,0.0631055,0.068826,0.0732845,0.076608,0.0789995,0.080673,0.081816,0.082576,0.0830655,0.083369,0.0835505,0.0836545,0.0837125,0.083744,0.08376,0.083769,0.083772,0.083774,0.083775,0.083776,0.083776,0.083776,0.083776,0.083776,
		0.015319,0.021761,0.030516,0.042078,0.056763,0.07447,0.094486,0.115519,0.135965,0.1543675,0.169799,0.181968,0.191081,0.197609,0.202097,0.2050535,0.2069165,0.2080385,0.2086865,0.2090475,0.209243,0.209344,0.209394,0.209419,0.209431,0.209436,0.209438,0.209439,0.209439,0.209439,0.209439,
		0.0163315,0.023575,0.033758,0.047842,0.0668895,0.0918635,0.1232365,0.1605515,0.2020905,0.2449705,0.285846,0.32188,0.3514215,0.374108,0.390509,0.4016975,0.4089075,0.413311,0.4158775,0.4173165,0.4180935,0.418496,0.418697,0.418796,0.418843,0.418865,0.418874,0.418877,0.418878,0.418879,0.418879,
		0.017074,0.0249075,0.036148,0.052124,0.074539,0.105451,0.1470675,0.2012445,0.268724,0.348099,0.4351855,0.523457,0.605704,0.6760965,0.7315175,0.771749,0.7988,0.8157805,0.8258375,0.83153,0.8346235,0.8362275,0.8370295,0.837424,0.837614,0.8377,0.837736,0.83775,0.837755,0.837757,0.837758,
		0.0177685,0.0261535,0.0383845,0.0561385,0.081742,0.118361,0.170169,0.242382,0.3410015,0.472022,0.639523,0.8427395,1.0729395,1.3122805,1.5370325,1.725589,1.8668195,1.962242,2.021453,2.055901,2.0749135,2.0848555,2.089847,2.092311,2.093497,2.094034,2.094259,2.094346,2.094378,2.094389,2.094393,
		0.018144,0.026828,0.039595,0.0583115,0.085646,0.1253765,0.18277,0.2650295,0.381667,0.5447375,0.768474,1.06739,1.451395,1.915988,2.430084,2.935145,3.367886,3.6915905,3.9068485,4.037592,4.1115845,4.1508055,4.1706365,4.180462,4.185199,4.187346,4.188244,4.188592,4.18872,4.188766,4.188782,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0111205,0.0145895,0.0185165,0.0226545,0.0266895,0.0303305,0.0333945,0.035825,0.0376635,0.039005,0.039957,0.0406185,0.0410705,0.0413735,0.041573,0.041701,0.04178,0.0418275,0.041855,0.0418705,0.041879,0.041883,0.041886,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.013067,0.01793,0.02403,0.03128,0.0393485,0.047676,0.0556185,0.062645,0.0684585,0.0730045,0.0764025,0.078853,0.080571,0.0817455,0.0825275,0.083032,0.083347,0.0835355,0.0836445,0.083706,0.08374,0.083758,0.083767,0.083772,0.083774,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,
		0.014895,0.0211845,0.0297505,0.0410925,0.0555465,0.073059,0.0929655,0.1139925,0.134542,0.153134,0.1687965,0.181195,0.1905085,0.1971985,0.201809,0.204855,0.20678,0.207946,0.208624,0.2090055,0.2092145,0.2093265,0.209384,0.209413,0.209428,0.209434,0.209438,0.209439,0.209439,0.209439,0.209439,
		0.01585,0.022902,0.032827,0.0465795,0.06522,0.0897345,0.1206655,0.15762,0.1989655,0.2418805,0.283017,0.3194685,0.349491,0.372639,0.3894355,0.4009345,0.408376,0.4129465,0.4156305,0.4171485,0.4179825,0.4184285,0.4186585,0.418774,0.418831,0.418859,0.418871,0.418876,0.418878,0.418879,0.418879,
		0.0165475,0.0241575,0.035085,0.050634,0.0724825,0.1026705,0.143424,0.196666,0.2632315,0.3418975,0.4286645,0.517111,0.599992,0.6713145,0.7277575,0.768938,0.7967745,0.81436,0.824863,0.8308645,0.834182,0.835958,0.8368765,0.8373385,0.837568,0.837677,0.837726,0.837746,0.837754,0.837757,0.837758,
		0.0171925,0.0253235,0.0371865,0.054416,0.0792845,0.11489,0.1653295,0.2357775,0.332221,0.460705,0.62555,0.8264205,1.05515,1.294386,1.520528,1.711597,1.8557915,1.9540465,2.015622,2.051836,2.0721875,2.083187,2.088893,2.0917745,2.093207,2.093891,2.094196,2.094321,2.094368,2.094386,2.094392,
		0.0175345,0.025946,0.0383155,0.0564525,0.0829535,0.121497,0.177221,0.2571685,0.370713,0.529742,0.7484085,1.041365,1.419046,1.8780875,2.3890875,2.8948915,3.332161,3.662689,3.885125,4.0219445,4.1008925,4.1441935,4.166839,4.1783225,4.18404,4.186775,4.187993,4.188492,4.188683,4.188752,4.188777,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.010933,0.0143745,0.0182795,0.0224095,0.026455,0.030126,0.0332285,0.035697,0.037569,0.0389375,0.039909,0.0405855,0.0410475,0.041358,0.041563,0.041694,0.0417755,0.0418245,0.041853,0.0418685,0.041878,0.041883,0.041885,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0128095,0.0176095,0.0236405,0.0308295,0.038859,0.047189,0.0551755,0.06227,0.068159,0.072777,0.076235,0.078734,0.0804875,0.0816875,0.082488,0.0830055,0.083329,0.083523,0.0836365,0.0837005,0.0837365,0.0837555,0.083766,0.083771,0.083773,0.083775,0.083775,0.083776,0.083776,0.083776,0.083776,
		0.0145615,0.020741,0.0291615,0.040332,0.0546015,0.071947,0.091756,0.112777,0.133406,0.152145,0.167989,0.18057,0.190045,0.196865,0.201575,0.2046925,0.206669,0.2078705,0.2085735,0.208972,0.209192,0.209312,0.2093755,0.209408,0.209425,0.209433,0.209437,0.209438,0.209439,0.209439,0.209439,
		0.015472,0.0223835,0.032114,0.0456105,0.0639335,0.088081,0.1186455,0.1553125,0.196495,0.239424,0.280755,0.317532,0.3479345,0.371451,0.3885645,0.4003145,0.407943,0.4126495,0.415429,0.4170135,0.417893,0.418369,0.4186225,0.418753,0.418819,0.418852,0.418868,0.418875,0.418877,0.418878,0.418879,
		0.0161325,0.0235785,0.03427,0.049493,0.070905,0.1005295,0.1405945,0.1930895,0.25893,0.3370135,0.423499,0.512057,0.5954195,0.6674695,0.724723,0.766661,0.795129,0.8132045,0.8240705,0.83033,0.8338225,0.8357245,0.836733,0.8372545,0.837519,0.837651,0.837713,0.837741,0.837752,0.837756,0.837757,
		0.0167365,0.024679,0.0362655,0.0530975,0.077404,0.11223,0.161612,0.23067,0.325409,0.4518975,0.614625,0.8135965,1.0410965,1.2801705,1.5073435,1.700358,1.8468915,1.947407,2.010896,2.04858,2.069969,2.0817305,2.088005,2.091253,2.0929035,2.093726,2.094115,2.094286,2.094354,2.094381,2.09439,
		0.017054,0.0252605,0.0373275,0.055026,0.080893,0.1185285,0.17297,0.2511335,0.362261,0.5181505,0.732849,1.0211075,1.393759,1.8483175,2.3567025,2.862889,3.303576,3.6394325,3.8676045,4.0094545,4.092218,4.1384375,4.1633095,4.17624,4.1828275,4.186113,4.18767,4.188352,4.188628,4.188732,4.18877,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.0107705,0.014206,0.018093,0.022215,0.0262655,0.0299555,0.033088,0.0355885,0.0374885,0.0388795,0.039869,0.040557,0.0410285,0.041345,0.0415535,0.041688,0.0417715,0.0418215,0.041851,0.0418675,0.041877,0.041882,0.041885,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0125895,0.017361,0.0233365,0.0304745,0.0384705,0.0467935,0.054809,0.061959,0.0679105,0.072588,0.0760965,0.0786345,0.080418,0.0816395,0.082455,0.0829835,0.0833135,0.083513,0.0836295,0.083696,0.083733,0.0837535,0.083764,0.08377,0.083773,0.083774,0.083775,0.083776,0.083776,0.083776,0.083776,
		0.0142835,0.020402,0.0287085,0.039742,0.053864,0.0710695,0.0907845,0.111796,0.132491,0.1513485,0.167338,0.1800655,0.1896705,0.196596,0.2013855,0.2045615,0.2065785,0.2078085,0.208532,0.208944,0.2091735,0.2092995,0.2093675,0.209403,0.209422,0.209431,0.209436,0.209438,0.209439,0.209439,0.209439,
		0.015165,0.0219945,0.0315725,0.044867,0.0629395,0.086796,0.117054,0.1534755,0.1945295,0.2374655,0.2789455,0.315977,0.3466815,0.370493,0.3878615,0.3998125,0.4075915,0.4124065,0.415264,0.416903,0.4178195,0.4183205,0.418591,0.4187335,0.418808,0.418845,0.418864,0.418873,0.418877,0.418878,0.418879,
		0.0158055,0.023155,0.0336635,0.0486315,0.069701,0.0988825,0.1384005,0.19028,0.25554,0.3331545,0.4193995,0.508029,0.591761,0.6643835,0.72228,0.7648215,0.793794,0.812263,0.823421,0.8298915,0.8335305,0.8355295,0.836607,0.8371775,0.8374725,0.837623,0.837698,0.837733,0.837748,0.837755,0.837757,
		0.016392,0.0242245,0.035601,0.052127,0.075999,0.110218,0.158773,0.2267365,0.320105,0.44502,0.606056,0.8034925,1.0299735,1.2688725,1.496817,1.691337,1.8397035,1.942012,2.007032,2.045911,2.068168,2.0805235,2.0872255,2.0907725,2.092612,2.093552,2.0940205,2.09424,2.094335,2.094373,2.094387,
		0.016691,0.024786,0.036631,0.053999,0.0793825,0.116319,0.1697665,0.2465415,0.3557665,0.5091695,0.7207515,1.0052875,1.373929,1.824875,2.331079,2.8374085,3.28065,3.620634,3.8533355,3.9992425,4.085182,4.1336715,4.1602105,4.174329,4.1816625,4.1854195,4.187292,4.188171,4.18855,4.188702,4.188758,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.010506,0.0140655,0.017948,0.022062,0.026114,0.0298155,0.032968,0.035494,0.037419,0.0388285,0.039833,0.0405325,0.0410115,0.041334,0.041546,0.0416825,0.041768,0.041819,0.041849,0.0418665,0.0418765,0.041882,0.041885,0.041886,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.012237,0.0171555,0.023101,0.030199,0.038166,0.046478,0.054508,0.0616975,0.0677005,0.072428,0.0759795,0.078551,0.08036,0.0815995,0.0824275,0.082964,0.0833005,0.0835045,0.0836235,0.0836915,0.0837305,0.0837515,0.083763,0.083769,0.083773,0.083774,0.083775,0.083775,0.083776,0.083776,0.083776,
		0.0138475,0.0201255,0.0283625,0.039289,0.053294,0.0703865,0.090014,0.111,0.1317455,0.150702,0.166811,0.179658,0.1893675,0.1963785,0.2012325,0.2044545,0.206505,0.2077585,0.2084975,0.208921,0.2091585,0.209289,0.209361,0.209399,0.209419,0.209429,0.209435,0.209437,0.209439,0.209439,0.209439,
		0.0146885,0.021682,0.031163,0.044303,0.0621795,0.0858055,0.115817,0.152018,0.192958,0.235904,0.2775025,0.314735,0.345679,0.369725,0.387296,0.399407,0.407306,0.41221,0.4151295,0.4168135,0.41776,0.418281,0.418565,0.418717,0.4187975,0.418839,0.41886,0.418871,0.418876,0.418878,0.418879,
		0.0153045,0.02282,0.0332115,0.0479845,0.06879,0.097626,0.1367145,0.188094,0.2528635,0.3301095,0.4161575,0.5048305,0.588848,0.6619195,0.7203215,0.7633395,0.792714,0.811499,0.8228945,0.8295365,0.833294,0.8353745,0.8365055,0.837112,0.837432,0.8375975,0.837682,0.837725,0.837744,0.837753,0.837756,
		0.0158695,0.023875,0.0351175,0.0514145,0.074955,0.1087075,0.156621,0.2237285,0.316,0.439643,0.599348,0.7955495,1.0211965,1.2599215,1.488432,1.684106,1.833912,1.937651,2.003908,2.043753,2.066717,2.0795605,2.0865905,2.0903675,2.0923575,2.093392,2.093923,2.094187,2.09431,2.094362,2.094383,
		0.0161545,0.024431,0.036138,0.053261,0.0782795,0.1146845,0.1673685,0.243067,0.350808,0.502228,0.711346,0.992948,1.358402,1.8064445,2.3108205,2.817122,3.2622815,3.6055095,3.841838,3.991004,4.0795215,4.12987,4.157688,4.1727155,4.180649,4.184778,4.1869015,4.187957,4.188449,4.188659,4.188742,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.009636,0.01381,0.017826,0.021947,0.026,0.0297065,0.0328695,0.035412,0.0373555,0.038783,0.0398005,0.0405105,0.0409965,0.041323,0.041539,0.041678,0.0417645,0.041817,0.0418475,0.0418655,0.0418755,0.0418815,0.0418845,0.041886,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.01111,0.016787,0.022905,0.029993,0.037939,0.0462395,0.0542715,0.0614815,0.067522,0.0722905,0.0758785,0.078479,0.0803095,0.0815645,0.0824035,0.0829475,0.0832895,0.0834965,0.0836185,0.0836885,0.083728,0.08375,0.0837625,0.083769,0.083772,0.083774,0.083775,0.083775,0.083776,0.083776,0.083776,
		0.012482,0.0196425,0.028076,0.038954,0.052874,0.0698785,0.0894305,0.110374,0.13114,0.1501725,0.1663795,0.1793245,0.1891205,0.1962,0.201106,0.2043665,0.206444,0.2077165,0.2084695,0.2089025,0.209146,0.2092805,0.209355,0.2093955,0.209417,0.209428,0.209434,0.209437,0.209438,0.209439,0.209439,
		0.0132015,0.0211375,0.0308225,0.0438845,0.0616215,0.085073,0.1148905,0.1509035,0.191716,0.2346595,0.276356,0.3137515,0.344885,0.3691155,0.3868445,0.3990815,0.407077,0.412051,0.415022,0.4167415,0.4177125,0.4182505,0.4185435,0.4187035,0.418789,0.4188335,0.418857,0.418869,0.418874,0.418877,0.418878,
		0.013719,0.022224,0.0328225,0.0474945,0.0681155,0.096696,0.1354555,0.186445,0.2507975,0.3277175,0.413617,0.5023245,0.5865575,0.659972,0.718764,0.762153,0.791846,0.810885,0.822472,0.8292515,0.833106,0.835251,0.8364255,0.8370585,0.8373985,0.8375765,0.8376685,0.837716,0.83774,0.837751,0.837755,
		0.0141645,0.0232065,0.034657,0.050829,0.074143,0.1075605,0.1549975,0.2214555,0.3128795,0.435488,0.5941295,0.789371,1.014343,1.2528915,1.481799,1.6783445,1.829276,1.9341565,2.0014045,2.042027,2.065557,2.0787945,2.0860885,2.0900425,2.0921495,2.093259,2.0938355,2.094132,2.09428,2.094348,2.094377,
		0.0143595,0.023687,0.0356035,0.0525885,0.077358,0.113386,0.1655135,0.2404095,0.3470255,0.4969135,0.704062,0.9833845,1.346341,1.792051,2.294881,2.801041,3.2476375,3.5934225,3.832646,3.98443,4.075004,4.126847,4.1556915,4.171414,4.179829,4.1842475,4.186551,4.187739,4.188329,4.188602,4.188718,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.007856,0.013075,0.0176255,0.0218475,0.0259125,0.0296235,0.0327915,0.0353425,0.037298,0.038739,0.039769,0.0404885,0.040981,0.0413125,0.041532,0.041673,0.041761,0.041815,0.041846,0.041864,0.041875,0.041881,0.041884,0.041886,0.041887,0.041887,0.041888,0.041888,0.041888,0.041888,0.041888,
		0.0088575,0.015744,0.022586,0.0298155,0.037767,0.046062,0.0540925,0.061311,0.067372,0.07217,0.0757885,0.0784145,0.080264,0.081533,0.082382,0.082933,0.083279,0.08349,0.083614,0.0836855,0.083726,0.0837485,0.0837615,0.0837685,0.083772,0.083774,0.083775,0.083775,0.083776,0.083776,0.083776,
		0.009756,0.018278,0.02762,0.0386675,0.052558,0.0695055,0.0890015,0.1099045,0.1306625,0.149737,0.1660195,0.179046,0.1889135,0.1960505,0.2009995,0.2042915,0.2063925,0.2076815,0.2084455,0.208887,0.209136,0.209274,0.2093505,0.2093925,0.2094155,0.209427,0.209433,0.209436,0.209438,0.209439,0.209439,
		0.010196,0.0195865,0.0302855,0.0435255,0.0612015,0.0845385,0.1142155,0.1500835,0.190774,0.2336735,0.275438,0.3129655,0.344252,0.368629,0.3864825,0.3988205,0.4068925,0.4119235,0.414936,0.416684,0.4176745,0.4182255,0.418527,0.418693,0.4187825,0.4188295,0.418854,0.418867,0.418873,0.418876,0.418878,
		0.010486,0.0205115,0.032213,0.0470695,0.0676005,0.0960145,0.1345415,0.18524,0.249267,0.325885,0.411637,0.5003805,0.584783,0.6584545,0.7175415,0.7612175,0.79116,0.8104005,0.8221395,0.829028,0.832958,0.8351545,0.836363,0.837017,0.837371,0.8375595,0.8376575,0.837708,0.837735,0.837748,0.837754,
		0.010718,0.021318,0.0339575,0.050317,0.073508,0.1067045,0.153811,0.2197975,0.310588,0.4323985,0.59015,0.7846375,1.009083,1.2474515,1.476621,1.6738165,1.8256195,1.9314015,1.999436,2.0406745,2.064652,2.0781965,2.085697,2.0897875,2.0919865,2.093156,2.0937675,2.094085,2.09425,2.094332,2.094369,
		0.0108275,0.0217245,0.034865,0.0520265,0.07664,0.112406,0.1641425,0.2384615,0.344248,0.4929865,0.6986105,0.9761,1.337129,1.780972,2.2824985,2.788458,3.2361365,3.5839235,3.825435,3.9792835,4.071481,4.124494,4.1541415,4.1703995,4.1791815,4.183838,4.1862795,4.18755,4.1882085,4.188536,4.188687};
	const double kb[]={							// 200 pts, eps=1e-6, down only
		// within a row steps along s in increasing order; going down rows steps along p in increasing order; going down groups steps along b in increasing order
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.04189,0.041893,0.041899,0.041915,0.041946,0.042004,0.042102,0.042253,0.042472,0.042791,0.043265,0.043982,0.045075,0.046771,0.049467,0.053936,0.061915,0.078412,0.127345,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083779,0.083783,0.083795,0.083822,0.083882,0.084007,0.084242,0.084639,0.085249,0.086139,0.087438,0.089387,0.092364,0.09699,0.104367,0.116654,0.138776,0.185554,0.336727,-1,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209441,0.209443,0.209448,0.209459,0.209487,0.209558,0.209727,0.210106,0.210891,0.212377,0.214913,0.218849,0.224647,0.233238,0.246406,0.267199,0.30123,0.360376,0.475543,0.76736,2.892618,-1,-1,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418881,0.418884,0.418893,0.418913,0.418958,0.41907,0.419352,0.420032,0.421554,0.424732,0.430793,0.441264,0.457809,0.482671,0.5204,0.580362,0.680613,0.860587,1.230851,2.299917,43.833248,-1,-1,-1,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.837761,0.837766,0.83778,0.837814,0.837893,0.838075,0.83852,0.839651,0.842386,0.848543,0.861495,0.886639,0.931319,1.004618,1.119604,1.303281,1.618464,2.215848,3.563967,8.579097,-1,-1,-1,-1,-1,
		2.094395,2.094395,2.094395,2.094396,2.094397,2.094401,2.094413,2.094445,2.09453,2.094743,2.095239,2.096378,2.099163,2.106266,2.123597,2.163102,2.248734,2.424674,2.768633,3.416744,4.622493,6.973167,12.431788,32.955733,-1,-1,-1,-1,-1,-1,-1,
		4.18879,4.188791,4.188791,4.188793,4.188798,4.188814,4.18886,4.188988,4.18933,4.190181,4.192165,4.196726,4.207891,4.236582,4.307234,4.472345,4.84993,5.714886,7.801409,13.611263,35.665157,249.781305,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.04189,0.041893,0.041899,0.041914,0.041945,0.042002,0.042098,0.042246,0.042462,0.042775,0.043241,0.043944,0.045013,0.046661,0.049244,0.05339,0.060278,0.072304,0.095008,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083779,0.083783,0.083794,0.083821,0.08388,0.084003,0.084232,0.084622,0.085221,0.086095,0.087372,0.089285,0.092202,0.096716,0.103855,0.115514,0.135526,0.172778,0.252913,0.489165,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209441,0.209443,0.209448,0.209459,0.209486,0.209556,0.209722,0.210094,0.210863,0.212319,0.214802,0.218661,0.224355,0.232792,0.245693,0.265997,0.299011,0.355514,0.461024,0.694048,1.476607,-1,-1,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418881,0.418884,0.418893,0.418912,0.418957,0.419066,0.419344,0.420011,0.421504,0.424618,0.430552,0.440799,0.457001,0.481379,0.518352,0.576902,0.674227,0.846849,1.19031,2.065794,7.467055,-1,-1,-1,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.837761,0.837766,0.837779,0.837813,0.837891,0.83807,0.838507,0.839618,0.842301,0.84834,0.861027,0.885623,0.929279,1.000875,1.113229,1.292387,1.59785,2.169833,3.42083,7.467558,-1,-1,-1,-1,-1,
		2.094395,2.094395,2.094395,2.094396,2.094397,2.094401,2.094412,2.094444,2.094528,2.094737,2.095224,2.096344,2.099081,2.106056,2.123063,2.161785,2.245562,2.417206,2.751409,3.378164,4.539069,6.787577,11.901473,29.639946,-1,-1,-1,-1,-1,-1,-1,
		4.18879,4.188791,4.188791,4.188793,4.188798,4.188814,4.188859,4.188985,4.189321,4.190158,4.192108,4.196591,4.20756,4.235731,4.305038,4.466731,4.835328,5.67457,7.675004,13.112382,32.696648,185.944927,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.04189,0.041892,0.041899,0.041913,0.041941,0.041995,0.042086,0.042225,0.04243,0.042727,0.043168,0.043831,0.044833,0.046359,0.048705,0.052354,0.058145,0.067631,0.083951,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083779,0.083783,0.083793,0.083819,0.083874,0.083989,0.084205,0.08457,0.085135,0.085963,0.087175,0.088983,0.091723,0.095928,0.102469,0.112838,0.129721,0.158358,0.210233,0.314892,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209441,0.209443,0.209447,0.209458,0.209484,0.209549,0.209707,0.210057,0.210779,0.212145,0.214472,0.218102,0.223484,0.231459,0.243578,0.262478,0.292734,0.342919,0.430524,0.595614,0.94847,1.922866,8.664989,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418881,0.418884,0.418892,0.41891,0.418953,0.419056,0.419319,0.419949,0.421355,0.424279,0.429837,0.439424,0.454609,0.477541,0.512272,0.566736,0.655845,0.80934,1.094723,1.695519,3.296618,11.995254,-1,-1,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.837761,0.837765,0.837778,0.83781,0.837884,0.838054,0.838468,0.839518,0.842051,0.847734,0.859637,0.88262,0.923268,0.989857,1.09445,1.260436,1.538499,2.042534,3.067459,5.652284,17.619775,-1,-1,-1,-1,
		2.094395,2.094395,2.094395,2.094396,2.094397,2.094401,2.094411,2.094441,2.094521,2.094719,2.095182,2.096243,2.098834,2.105427,2.121469,2.157867,2.236183,2.395286,2.701418,3.267611,4.302655,6.270157,10.509836,22.764819,151.045632,-1,-1,-1,-1,-1,-1,
		4.18879,4.188791,4.188791,4.188793,4.188798,4.188813,4.188855,4.188975,4.189294,4.190087,4.191938,4.196186,4.20657,4.233189,4.298489,4.450081,4.792394,5.558028,7.321584,11.807962,25.956103,99.819378,-1,-1,-1,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.04189,0.041892,0.041898,0.04191,0.041936,0.041984,0.042065,0.04219,0.042376,0.042647,0.043048,0.043646,0.044545,0.045904,0.047975,0.051171,0.056188,0.064266,0.077746,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083778,0.083782,0.083792,0.083814,0.083865,0.083968,0.08416,0.084487,0.084994,0.085746,0.086848,0.088485,0.090946,0.094689,0.100439,0.109411,0.123741,0.14742,0.188511,0.265138,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209441,0.209442,0.209447,0.209456,0.20948,0.209539,0.209681,0.209996,0.210643,0.21186,0.213936,0.217188,0.222048,0.22926,0.24013,0.256861,0.283168,0.325616,0.39671,0.52261,0.765195,1.300398,2.843005,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418881,0.418884,0.418891,0.418908,0.418946,0.41904,0.419278,0.419846,0.421109,0.423725,0.428674,0.437192,0.450721,0.471268,0.502351,0.550442,0.627378,0.755435,0.979878,1.404128,2.30864,4.729251,17.562992,-1,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.83776,0.837765,0.837776,0.837805,0.837873,0.838027,0.838403,0.839353,0.841636,0.846741,0.85737,0.877755,0.913598,0.972195,1.064312,1.209536,1.446946,1.858685,2.632324,4.265466,8.499595,27.659554,-1,-1,-1,
		2.094395,2.094395,2.094395,2.094396,2.094397,2.0944,2.09441,2.094437,2.09451,2.09469,2.095111,2.096075,2.098425,2.104388,2.118842,2.151452,2.220979,2.360289,2.623308,3.099209,3.950615,5.522913,8.685969,16.262888,43.774442,-1,-1,-1,-1,-1,-1,
		4.18879,4.188791,4.188791,4.188793,4.188797,4.188811,4.18885,4.188958,4.189248,4.18997,4.191654,4.195514,4.204928,4.228987,4.287729,4.422951,4.723633,5.377343,6.806611,10.118009,18.968705,51.250113,463.883025,-1,-1,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.041891,0.041896,0.041907,0.041928,0.041969,0.042037,0.042143,0.042302,0.042536,0.042881,0.043393,0.044158,0.04531,0.047055,0.049723,0.053854,0.060368,0.07091,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083778,0.083781,0.083789,0.083809,0.083851,0.083938,0.0841,0.084373,0.084802,0.085445,0.086394,0.087797,0.08989,0.093046,0.09785,0.105255,0.116871,0.135538,0.166612,0.221056,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209442,0.209446,0.209454,0.209474,0.209525,0.209646,0.209912,0.210457,0.211475,0.21321,0.215946,0.220078,0.226231,0.235437,0.249402,0.270965,0.304972,0.360171,0.453538,0.621092,0.949585,1.693391,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.41888,0.418883,0.418889,0.418904,0.418937,0.419017,0.41922,0.419704,0.420773,0.422971,0.427104,0.434192,0.445488,0.462765,0.488911,0.528818,0.591072,0.691204,0.858741,1.15507,1.724105,2.967727,6.44628,28.575918,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.83776,0.837764,0.837774,0.837799,0.837857,0.837989,0.838312,0.839124,0.841065,0.845381,0.854293,0.871226,0.900757,0.94888,1.024552,1.143004,1.331737,1.643974,2.190723,3.226575,5.438865,11.324272,40.751615,-1,-1,
		2.094395,2.094395,2.094395,2.094396,2.094397,2.0944,2.094408,2.094431,2.094494,2.094649,2.095012,2.095842,2.097857,2.102947,2.115228,2.142705,2.200544,2.314289,2.523772,2.892267,3.532851,4.673434,6.817306,11.253652,22.251893,63.769736,-1,-1,-1,-1,-1,
		4.18879,4.188791,4.188791,4.188792,4.188796,4.188808,4.188841,4.188935,4.189185,4.189807,4.191257,4.194578,4.202646,4.223175,4.272963,4.386229,4.632847,5.149415,6.208689,8.419692,13.446855,27.363129,86.948411,-1,-1,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.041891,0.041894,0.041903,0.041919,0.04195,0.042002,0.042084,0.042209,0.042394,0.042668,0.043073,0.043675,0.044575,0.04593,0.047979,0.051102,0.055915,0.063458,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083778,0.08378,0.083787,0.083802,0.083835,0.083901,0.084025,0.084234,0.084564,0.085067,0.085816,0.086926,0.08857,0.091026,0.094726,0.10035,0.108989,0.122466,0.143985,0.179561,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209442,0.209444,0.209451,0.209467,0.209507,0.209601,0.209808,0.210226,0.211002,0.212322,0.214418,0.217624,0.222435,0.229608,0.240345,0.256614,0.281666,0.32101,0.38452,0.491294,0.682052,1.055726,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418882,0.418887,0.418899,0.418925,0.418989,0.419148,0.419525,0.420352,0.422038,0.425179,0.430542,0.439118,0.452348,0.47243,0.502789,0.549049,0.620952,0.736057,0.927759,1.265585,1.913555,3.335807,7.380686,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.83776,0.837763,0.837771,0.83779,0.837836,0.837942,0.838196,0.838834,0.840346,0.843683,0.850498,0.863286,0.885362,0.92119,0.977526,1.065184,1.201858,1.418619,1.775159,2.39651,3.572477,6.094668,12.893471,48.988662,-1,
		2.094395,2.094395,2.094395,2.094396,2.094396,2.094399,2.094405,2.094424,2.094473,2.094597,2.094885,2.095543,2.097133,2.101122,2.110685,2.131845,2.175647,2.259799,2.410317,2.666876,3.098671,3.837355,5.144837,7.589286,12.616051,25.008685,71.833542,-1,-1,-1,-1,
		4.18879,4.188791,4.188791,4.188792,4.188795,4.188804,4.188831,4.188905,4.189103,4.189598,4.190749,4.19338,4.199742,4.215829,4.254483,4.341072,4.524652,4.892536,5.596009,6.923069,9.549078,15.456381,32.17364,107.399061,-1,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.04189,0.041892,0.041898,0.041909,0.041928,0.041962,0.042015,0.042098,0.042224,0.042411,0.042688,0.043098,0.043706,0.044615,0.045974,0.048014,0.051094,0.055783,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083779,0.083783,0.083793,0.083815,0.083858,0.083938,0.084072,0.084287,0.08462,0.085126,0.085881,0.086999,0.088655,0.091121,0.094814,0.100377,0.108826,0.121827,0.142276,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209441,0.209443,0.209447,0.209458,0.209485,0.209548,0.209684,0.209956,0.210454,0.2113,0.212655,0.214763,0.217972,0.222781,0.229934,0.240594,0.256616,0.281005,0.318774,0.378813,0.478195,0.652985,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418881,0.418885,0.418892,0.41891,0.418954,0.419061,0.419312,0.419857,0.420948,0.42296,0.426375,0.431866,0.440443,0.453585,0.473442,0.503276,0.54834,0.617554,0.726783,0.905991,1.217112,1.804415,3.066873,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.837761,0.837767,0.83778,0.837812,0.837884,0.838056,0.838485,0.83949,0.841682,0.84609,0.854224,0.8681,0.890565,0.926018,0.981182,1.066264,1.19756,1.403302,1.737199,2.31088,3.380123,5.629203,11.488793,38.919817,
		2.094395,2.094395,2.094395,2.094395,2.094396,2.094397,2.094402,2.094415,2.094449,2.094533,2.09473,2.09518,2.096258,2.098935,2.10529,2.119132,2.147154,2.199476,2.290164,2.440395,2.686682,3.094438,3.783586,4.988204,7.210043,11.701616,22.465684,60.214519,-1,-1,-1,
		4.18879,4.18879,4.188791,4.188791,4.188793,4.1888,4.188818,4.188869,4.189004,4.189343,4.190131,4.191927,4.196235,4.207031,4.23263,4.288803,4.403991,4.62352,5.015335,5.695149,6.911447,9.276609,14.513104,28.774541,86.608354,-1,-1,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041889,0.04189,0.041892,0.041896,0.041904,0.041917,0.041939,0.041974,0.042029,0.042113,0.042241,0.042431,0.042712,0.04313,0.043756,0.044695,0.046111,0.048261,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083777,0.083779,0.083783,0.083792,0.08381,0.083841,0.083894,0.08398,0.08412,0.084341,0.084683,0.085198,0.085967,0.08711,0.088817,0.091381,0.095265,0.101197,0.110371,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209441,0.209443,0.209448,0.20946,0.209487,0.209544,0.209653,0.209849,0.210178,0.21072,0.211601,0.213,0.215169,0.218466,0.223421,0.230854,0.242079,0.259246,0.28593,0.328323,0.397876,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.41888,0.418882,0.418885,0.418893,0.418913,0.41896,0.419068,0.419294,0.419732,0.420516,0.42184,0.424022,0.427577,0.433263,0.442155,0.45583,0.476659,0.508391,0.557314,0.634483,0.760212,0.974604,1.3659,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837758,0.837759,0.837759,0.837762,0.837768,0.837783,0.837816,0.837893,0.838082,0.838511,0.839418,0.841186,0.844347,0.849676,0.858474,0.872933,0.89639,0.933811,0.992776,1.0852,1.23095,1.466267,1.86363,2.581515,4.0141,7.393281,
		2.094395,2.094395,2.094395,2.094395,2.094395,2.094396,2.094398,2.094404,2.094419,2.094458,2.094549,2.094754,2.095238,2.096407,2.099126,2.104858,2.115981,2.135879,2.169474,2.225581,2.320553,2.482271,2.757847,3.230777,4.057986,5.558948,8.47357,14.919604,33.625361,178.666091,-1,
		4.18879,4.18879,4.18879,4.188791,4.188792,4.188794,4.188803,4.188826,4.188887,4.189042,4.189404,4.190225,4.192152,4.196886,4.207805,4.230858,4.27568,4.355961,4.49217,4.725153,5.14315,5.92393,7.44638,10.669479,18.625516,44.765654,235.432933,-1,-1,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041887,0.041886,0.041883,0.041878,0.04187,0.041857,0.04184,0.041817,0.041786,0.041744,0.041688,0.041614,0.041517,0.0414,0.041277,0.041192,0.041357,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083775,0.083775,0.083772,0.083768,0.083757,0.083738,0.083704,0.083654,0.083585,0.083494,0.083375,0.083223,0.083034,0.082806,0.082562,0.082365,0.082359,0.082834,0.084164,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.20944,0.209439,0.209439,0.209439,0.209439,0.209439,0.209438,0.209436,0.209431,0.209419,0.209389,0.209325,0.209201,0.208993,0.208683,0.208267,0.207742,0.207114,0.20642,0.205763,0.205369,0.205635,0.207466,0.212094,0.221526,0.238801,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418878,0.418877,0.418874,0.418867,0.418846,0.418795,0.418676,0.418418,0.417923,0.417096,0.415896,0.414343,0.412525,0.410667,0.40925,0.409173,0.411945,0.419919,0.436626,0.467499,0.521432,0.61418,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837758,0.837758,0.837757,0.837757,0.837754,0.837749,0.837738,0.837708,0.837627,0.837421,0.836937,0.835908,0.833959,0.830756,0.826278,0.820945,0.815734,0.81252,0.814517,0.826689,0.856202,0.913374,1.014188,1.186764,1.486929,2.037585,
		2.094395,2.094395,2.094395,2.094395,2.094395,2.094395,2.094394,2.094392,2.094386,2.094372,2.09434,2.094268,2.094081,2.093565,2.092287,2.089324,2.083029,2.071286,2.052973,2.030307,2.010196,2.004394,2.03042,2.113761,2.293162,2.633395,3.25431,4.401342,6.652126,11.715432,27.281709,
		4.18879,4.18879,4.18879,4.18879,4.188789,4.188788,4.188785,4.188776,4.188753,4.188697,4.188571,4.18828,4.18752,4.185514,4.180434,4.168698,4.144066,4.099671,4.036144,3.973002,3.950515,4.023319,4.267427,4.815758,5.953117,8.408489,14.49376,34.595397,198.536869,-1,-1,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041887,0.041886,0.041884,0.041879,0.041869,0.041851,0.041821,0.041773,0.041701,0.041593,0.041437,0.041209,0.040883,0.040422,0.039783,0.038918,0.037793,0.036398,0.034773,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083775,0.083775,0.083773,0.083769,0.083761,0.083741,0.083702,0.083629,0.083508,0.083318,0.083031,0.082611,0.082004,0.081143,0.079943,0.078315,0.07618,0.073507,0.070347,0.066865,0.063332,
		0.20944,0.20944,0.20944,0.20944,0.20944,0.209439,0.209439,0.209439,0.209439,0.209438,0.209437,0.209433,0.209424,0.2094,0.209344,0.209223,0.208977,0.208527,0.207776,0.20661,0.204886,0.202421,0.198999,0.194403,0.188475,0.181209,0.172846,0.163915,0.155166,0.147438,0.141538,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418878,0.418877,0.418874,0.418867,0.418852,0.418815,0.418721,0.418498,0.418011,0.417032,0.415246,0.412294,0.407793,0.401324,0.392445,0.380813,0.36636,0.349494,0.331216,0.313027,0.296664,0.283811,0.276014,0.274684,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837758,0.837757,0.837755,0.83775,0.837739,0.837712,0.837651,0.837502,0.837126,0.836234,0.834288,0.830393,0.82336,0.811965,0.795204,0.772331,0.743144,0.708384,0.670026,0.631146,0.595394,0.566404,0.547468,0.541857,0.554099,0.590472,
		2.094395,2.094395,2.094395,2.094395,2.094394,2.094393,2.094389,2.094378,2.094348,2.094275,2.094106,2.093723,2.092792,2.09044,2.084874,2.072842,2.049137,2.007585,1.943899,1.858302,1.75572,1.644561,1.535149,1.437792,1.361711,1.315599,1.30964,1.358984,1.48978,1.750887,2.24019,
		4.18879,4.18879,4.18879,4.188789,4.188787,4.188781,4.188765,4.18872,4.188601,4.188309,4.187635,4.186101,4.18237,4.173031,4.150961,4.10373,4.012828,3.860917,3.647052,3.39409,3.135636,2.899599,2.705441,2.568647,2.506425,2.545535,2.73689,3.18998,4.170469,6.435561,12.763484,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041888,0.041887,0.041887,0.041885,0.041881,0.041872,0.041855,0.041824,0.041771,0.041688,0.041561,0.04137,0.041086,0.040671,0.040071,0.039222,0.038046,0.036469,0.034439,0.031955,0.029089,
		0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083776,0.083775,0.083775,0.083774,0.083771,0.083764,0.083748,0.083713,0.083644,0.083519,0.08331,0.082981,0.08248,0.081734,0.080646,0.07909,0.076917,0.073973,0.070138,0.065377,0.059797,0.05365,0.047287,
		0.20944,0.20944,0.20944,0.20944,0.209439,0.209439,0.209439,0.209439,0.209439,0.209437,0.209434,0.209427,0.209409,0.209366,0.209265,0.209047,0.208617,0.207842,0.206555,0.204554,0.20157,0.197269,0.191265,0.183187,0.172789,0.160092,0.145489,0.129723,0.113701,0.098262,0.083997,
		0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418879,0.418878,0.418875,0.41887,0.418857,0.418828,0.418758,0.418585,0.41818,0.417312,0.415603,0.412538,0.407522,0.399891,0.388889,0.373749,0.353906,0.329279,0.300516,0.269002,0.236562,0.204959,0.175514,0.148984,0.125658,
		0.837758,0.837758,0.837758,0.837758,0.837758,0.837757,0.837756,0.837753,0.837744,0.837722,0.83767,0.837554,0.837276,0.836584,0.834966,0.831508,0.82475,0.812819,0.793796,0.766027,0.728291,0.68028,0.623153,0.559645,0.493504,0.428482,0.367482,0.312235,0.263466,0.221213,0.185128,
		2.094395,2.094395,2.094395,2.094395,2.094394,2.094391,2.094383,2.094362,2.094306,2.094167,2.093847,2.093121,2.091381,2.08706,2.076991,2.055714,2.015059,1.946236,1.844179,1.710273,1.551713,1.379401,1.205005,1.038065,0.884841,0.748227,0.628948,0.526405,0.439294,0.366007,0.30487,
		4.18879,4.18879,4.18879,4.188788,4.188784,4.188774,4.188743,4.188658,4.188433,4.187877,4.186596,4.183696,4.176735,4.159569,4.119833,4.037265,3.884886,3.643077,3.319535,2.948997,2.570019,2.208704,1.878768,1.585867,1.330887,1.112094,0.926487,0.77048,0.640347,0.532576,0.444031};
	const double slo=-3,sinc=0.2,shi=3;						// logs of values; shi=3
	const double blor=0,bincr=0.2;								// logs of values; bhir=3
	const double blob=0,bincb=0.1;								// actual values; bhib=1
	const double pvect[]={0,0.01,0.02,0.05,0.1,0.2,0.5,1.0};
	const int snum=31,bnumr=16,bnumb=11,pnum=8;

	int sindx,bindx,pindx,sflag;
	double logstep,steplo,stephi,problo,probhi,blo,bhi,step0,step1,logb;
	double denom,numll,numlh,numhl,numhh,kll,klh,khl,khh;
	double numlll,numllh,numlhl,numlhh,numhll,numhlh,numhhl,numhhh;
	double klll,kllh,klhl,klhh,khll,khlh,khhl,khhh;
	double lambdap,theory,ans,ansp1;

	if(prob==1) return numrxnrate(step,a,b);					// reaction probability of 1 is better handled by prior function
	if(step<0 || a<0) return -1;											// rms step length and binding radius need to be >0
	if(prob<0 || prob>1) return -1;										// reaction probability cannot be outside of [0,1]
	if(prob==0) return 0;															// reaction probability of 0 has 0 reaction rate
	if(a==0) return 0;																// binding radius of 0 has 0 reaction rate
	if(step==0 && b>=0 && b<a) return -1;							// reaction rate is infinite with small steps and unbinding radius inside binding radius
	if(step==0) return 0;															// reaction rate is 0 if rms step length is 0
	// now, inputs are 0 < step < infinity, 0 < a < infinity, -infinity < b < infinity, 0 < prob < 1

	ansp1=numrxnrate(step,a,b);
	step/=a;																					// step is now reduced step length
	logstep=log(step);																// log of reduced step length
	b/=a;																							// reduced unbinding radius (or negative for irreversible)

	if(b<0 || b>=1.0) {
		sindx=(int)floor((logstep-slo)/sinc);						// find step index value, constrained by table
		if(sindx<0) {sindx=0;sflag=-1; }
		else if(sindx>snum-2) sindx=snum-2;
		else sflag=0;
		steplo=slo+sindx*sinc;													// tabulated step values low and high
		stephi=slo+(sindx+1)*sinc; }
	else {
		sindx=(int)floor((shi-logstep)/sinc);						// this table uses reversed step values
		if(sindx<0) {sindx=0;sflag=-1; }
		else if(sindx>snum-2) sindx=snum-2;
		else sflag=0;
		steplo=shi-sindx*sinc;													// tabulated step values low and high
		stephi=shi-(sindx+1)*sinc; }

	for(pindx=0;pindx<pnum-2 && prob>pvect[pindx+1];pindx++);
	problo=pvect[pindx];
	probhi=pvect[pindx+1];

	if(b<0) {																					// irreversible
		denom=(probhi-problo)*(stephi-steplo);
		numll=(probhi-prob)*(stephi-logstep);
		numlh=(probhi-prob)*(logstep-steplo);
		numhl=(prob-problo)*(stephi-logstep);
		numhh=(prob-problo)*(logstep-steplo);
		kll=kirr[pindx*snum+sindx];
		klh=kirr[pindx*snum+sindx+1];
		khl=kirr[(pindx+1)*snum+sindx];
		khh=kirr[(pindx+1)*snum+sindx+1];
		ans=(numll*kll+numlh*klh+numhl*khl+numhh*khh)/denom;
		if(sindx==0) {
			lambdap=sqrt(2*prob)/step;
			theory=2*PI*step*step*(1-tanh(lambdap)/lambdap);		// from Notes18C06 lambda-rho algorithm
			if(sflag==-1) ans=theory;														// use theory if step is below table
			else {																							// interpolate between theory and table if in lowest region
				step0=exp(steplo);
				step1=exp(steplo+sinc);
				ans=((step-step0)*ans+(step1-step)*theory)/(step1-step0); }}
		else if(ans>4.0*PI/3.0*prob)
			ans=4.0*PI/3.0*prob; }

	else if(b>=1.0) {																				// unbinding radius >= binding radius
		logb=log(b);
		bindx=(int)floor((logb-blor)/bincr);
		if(bindx<0) bindx=0;
		else if(bindx>bnumr-2) bindx=bnumr-2;
		blo=blor+bindx*bincr;
		bhi=blor+(bindx+1)*bincr;
		denom=(bhi-blo)*(probhi-problo)*(stephi-steplo);
		numlll=(bhi-logb)*(probhi-prob)*(stephi-logstep);
		numllh=(bhi-logb)*(probhi-prob)*(logstep-steplo);
		numlhl=(bhi-logb)*(prob-problo)*(stephi-logstep);
		numlhh=(bhi-logb)*(prob-problo)*(logstep-steplo);
		numhll=(logb-blo)*(probhi-prob)*(stephi-logstep);
		numhlh=(logb-blo)*(probhi-prob)*(logstep-steplo);
		numhhl=(logb-blo)*(prob-problo)*(stephi-logstep);
		numhhh=(logb-blo)*(prob-problo)*(logstep-steplo);
		klll=krev[bindx*snum*pnum+pindx*snum+sindx];
		kllh=krev[bindx*snum*pnum+pindx*snum+sindx+1];
		klhl=krev[bindx*snum*pnum+(pindx+1)*snum+sindx];
		klhh=krev[bindx*snum*pnum+(pindx+1)*snum+sindx+1];
		khll=krev[(bindx+1)*snum*pnum+pindx*snum+sindx];
		khlh=krev[(bindx+1)*snum*pnum+pindx*snum+sindx+1];
		khhl=krev[(bindx+1)*snum*pnum+(pindx+1)*snum+sindx];
		khhh=krev[(bindx+1)*snum*pnum+(pindx+1)*snum+sindx+1];
		ans=(numlll*klll+numllh*kllh+numlhl*klhl+numlhh*klhh+numhll*khll+numhlh*khlh+numhhl*khhl+numhhh*khhh)/denom;
		if(sindx==0) {
			lambdap=sqrt(2*prob)/step;
			if(lambdap<10) theory=2*PI*step*step*b*(lambdap*cosh(lambdap)-sinh(lambdap))/(lambdap*(b-1)*cosh(lambdap)+sinh(lambdap));
			else theory=2*PI*step*step*b*(lambdap-1)/(lambdap*(b-1)+1);
			if(sflag==-1) ans=theory;														// use theory if step is below table
			else {																							// interpolate between theory and table if in lowest region
				step0=exp(steplo);
				step1=exp(steplo+sinc);
				ans=((step-step0)*ans+(step1-step)*theory)/(step1-step0); }}
		else if(ans>4.0*PI/3.0*prob)
			ans=4.0*PI/3.0*prob; }

	else {																									// unbinding radius < binding radius
		bindx=(int)floor((b-blob)/bincb);
		if(bindx<0) bindx=0;
		else if(bindx>bnumb-2) bindx=bnumb-2;
		blo=blob+bindx*bincb;
		bhi=blob+(bindx+1)*bincb;
		while(kb[bindx*snum*pnum+(pindx+1)*snum+sindx+1]==-1) sindx--;		// avoid bad data values
		steplo=shi-sindx*sinc;													// tabulated step values low and high
		stephi=shi-(sindx+1)*sinc;
		denom=(bhi-blo)*(probhi-problo)*(stephi-steplo);
		numlll=(bhi-b)*(probhi-prob)*(stephi-logstep);
		numllh=(bhi-b)*(probhi-prob)*(logstep-steplo);
		numlhl=(bhi-b)*(prob-problo)*(stephi-logstep);
		numlhh=(bhi-b)*(prob-problo)*(logstep-steplo);
		numhll=(b-blo)*(probhi-prob)*(stephi-logstep);
		numhlh=(b-blo)*(probhi-prob)*(logstep-steplo);
		numhhl=(b-blo)*(prob-problo)*(stephi-logstep);
		numhhh=(b-blo)*(prob-problo)*(logstep-steplo);
		klll=kb[bindx*snum*pnum+pindx*snum+sindx];
		kllh=kb[bindx*snum*pnum+pindx*snum+sindx+1];
		klhl=kb[bindx*snum*pnum+(pindx+1)*snum+sindx];
		klhh=kb[bindx*snum*pnum+(pindx+1)*snum+sindx+1];
		khll=kb[(bindx+1)*snum*pnum+pindx*snum+sindx];
		khlh=kb[(bindx+1)*snum*pnum+pindx*snum+sindx+1];
		khhl=kb[(bindx+1)*snum*pnum+(pindx+1)*snum+sindx];
		khhh=kb[(bindx+1)*snum*pnum+(pindx+1)*snum+sindx+1];
		ans=(numlll*klll+numllh*kllh+numlhl*klhl+numlhh*klhh+numhll*khll+numhlh*khlh+numhhl*khhl+numhhh*khhh)/denom; }

	ans*=a*a*a;
	if(ans>ansp1) ans=ansp1;				// limit maximum value to that if probability is 1

	return ans; }


/* actrxnrate calculates the effective activation limited reaction rate for the
simulation, which is the reaction rate if the radial correlation function is 1
for all r>a.  The returned value needs to be divided by delta_t.  The equation
is ka=4*pi/3*[erfc(sqrt(2)/s)+s*sqrt((2/pi))]+2*sqrt(2*pi)/3*s*(s^2-1)*[exp(-2/s^2)-1].
It was calculated analytically and verified numerically. */
double actrxnrate(double step,double a,double prob) {
	double ka;

	if(step<0 || a<=0) return -1;
	if(step==0) return 0;
	step/=a;
	ka=4.0*PI/3.0*(rxnparam_erfccD(sqrt(2.0)/step)+step*sqrt(2.0/PI));
	ka+=2.0*SQRT2PI/3.0*step*(step*step-1.0)*(exp(-2.0/step/step)-1.0);
	return prob*ka*a*a*a; }


/* bindingradius returns the binding radius that corresponds to some given
information.  rate is the actual rate constant (not reduced), dt is the time
step, and difc is the mutual diffusion constant (sum of reactant diffusion
constants).  If b is -1, the reaction is assumed to be irreversible; if b>=0 and
rel=0, then the b value is used as the unbinding radius; and if b>=0 and rel=1,
then the b value is used as the ratio of the unbinding to binding radius, b/a.
This algorithm executes a simple search from numrxnrate, based on the fact that
reaction rates monotonically increase with increasing a, for all the b value
possibilities.  The return value is usually the binding radius.  However, a
value of -1 signifies illegal input parameters.
Modified 2/22/08 to allow for dt==0. */
double bindingradius(double rate,double dt,double difc,double b,int rel) {
	double a,lo,dif,step;
	int n;

	if(rate<0 || dt<0 || difc<=0) return -1;
	if(rate==0) return 0;
	if(dt==0) {
		if(b<0) return rate/(4*PI*difc);
		if(rel && b>1) return rate*(1-1/b)/(4*PI*difc);
		if(rel && b<=1) return -1;
		if(b>0) return rate/(4*PI*difc+rate/b);
		return -1; }

	step=sqrt(2.0*difc*dt);
	lo=0;
	a=step;
	while(numrxnrate(step,a,rel?b*a:b)<rate*dt) {
		lo=a;
		a*=2.0; }
	dif=a-lo;
	for(n=0;n<15;n++) {
		dif*=0.5;
		a=lo+dif;
		if(numrxnrate(step,a,rel?b*a:b)<rate*dt) lo=a; }
	a=lo+0.5*dif;
	return a; }


/* bindingradiusprob */
double bindingradiusprob(double rate,double dt,double difc,double b,int rel,double chi,double *probptr) {
	double a,prob,step,lambdap,lambda;
	double xtry[2],xmult[2],xlo[2],xhi[2],errtry,errold,chitry,ratetry,errbest,abest,pbest,unbindrad;
	int it,indx;

	if(rate<0 || dt<0 || difc<=0 || chi>1) return -1;
	if(chi<0 && (!probptr || *probptr==1)) return bindingradius(rate,dt,difc,b,rel);
	if(chi<0 && probptr && *probptr==0) return 0;
	if(rate==0) {
		if(probptr && *probptr<0) *probptr=0;
		return 0; }
	// By now, rate>0, dt>=0, difc>0, b anything, rel 0 or 1, chi >=0 or *probptr in (0,1), probptr various

	if(dt==0) {
		if(chi<0) return -1;
		else if(b<0) a=rate/(4*PI*difc*chi);									// k = 4πD sigma_b chi  from general steady-state soln
		else if(rel && b>1) a=rate*(1-1/b)/(4*PI*difc*chi);
		else if(rel && b<=1) return -1;
		else if(b>0) a=rate*b/(rate+4*PI*difc*b*chi);						// k = 4πD sigma_b sigma_u chi / (sigma_u - sigma_b)  from general steady-state soln
		else return -1;
		lambdap=(3.18243*chi-3.40864*chi*chi+0.984385*chi*chi*chi)/pow(1-chi,2.08761);
		lambda=difc*lambdap/a*lambdap/a;
		if(probptr) *probptr=lambda;
		return a; }

	if(chi<0) {																							// unknown chi but known probability
		prob=pbest=*probptr;
		step=sqrt(2*difc*dt);
		xtry[0]=abest=step;																		// binding radius is being searched
		xmult[0]=0.1;
		xlo[0]=1.0e-15;
		xhi[0]=1.0e15;
		unbindrad=(b>=0)?((rel)?b*xtry[0]:b):-1;
		ratetry=numrxnrateprob(step,xtry[0],unbindrad,prob)/dt;
		errtry=errbest=(ratetry/rate-1.0)*(ratetry/rate-1.0);
		for(it=0;it<200 && errtry>1.0e-6;it++) {
			errold=errtry;
//			printf("a=%g rate=%g error=%g mult=%g \n",xtry[0],ratetry,errtry,xmult[0]); //?? DEBUG CODE
			xtry[0]=xtry[0]*(1.0+xmult[0]);
			if(xtry[0]<xlo[0]) xtry[0]=xlo[0];
			else if(xtry[0]>xhi[0]) xtry[0]=xhi[0];
			unbindrad=(b>=0)?((rel)?b*xtry[0]:b):-1;
			ratetry=numrxnrateprob(step,xtry[0],unbindrad,prob)/dt;
			errtry=(ratetry/rate-1.0)*(ratetry/rate-1.0);
			if(errtry<errold) 																	// better so go faster in same direction
				xmult[0]*=1.1;
			else																								// worse so slow down and reverse
				xmult[0]*=-0.5;
			if(errtry<errbest) {
				errbest=errtry;
				abest=xtry[0]; }}}

	else {
		step=sqrt(2*difc*dt);
		xtry[0]=abest=bindingradius(rate,dt,difc,b,rel);			// first try with p=1 and initialize
		xtry[1]=pbest=1;
		ratetry=rate;
		unbindrad=(b>=0)?((rel)?b*xtry[0]:b):-1;
		chitry=(unbindrad>=0) ? rate/(4*PI*difc*xtry[0]*unbindrad/(unbindrad-xtry[0])) : rate/(4*PI*difc*xtry[0]);
		if(chitry>chi) {																			// it's possible to improve
			xmult[0]=xmult[1]=0.1;
			xlo[0]=xlo[1]=1.0e-15;
			xhi[0]=1.0e15;
			xhi[1]=1.0;
			xtry[0]=abest=step;
			xtry[1]=pbest=0.5;
			ratetry=numrxnrateprob(step,xtry[0],rel?b*xtry[0]:b,xtry[1])/dt;
			unbindrad=(b>=0)?((rel)?b*xtry[0]:b):-1;
			chitry=(unbindrad>=0) ? rate/(4*PI*difc*xtry[0]*unbindrad/(unbindrad-xtry[0])) : rate/(4*PI*difc*xtry[0]);
			errtry=errbest=(ratetry/rate-1.0)*(ratetry/rate-1.0)+(chitry/chi-1.0)*(chitry/chi-1.0);
			for(it=0;it<300 && errtry>1.0e-6;it++) {
				for(indx=0;indx<2;indx++) {
					errold=errtry;
//					printf("a %g prob %g rate %g chi %g error %g a_mult %g p_mult %g\n",xtry[0],xtry[1],ratetry,chitry,errtry,xmult[0],xmult[1]); //?? DEBUG CODE
					xtry[indx]=xtry[indx]*(1.0+xmult[indx]);
					if(xtry[indx]<xlo[indx]) xtry[indx]=xlo[indx];
					else if(xtry[indx]>xhi[indx]) xtry[indx]=xhi[indx];
					ratetry=numrxnrateprob(step,xtry[0],rel?b*xtry[0]:b,xtry[1])/dt;
					unbindrad=(b>=0)?((rel)?b*xtry[0]:b):-1;
					chitry=(unbindrad>=0) ? rate/(4*PI*difc*xtry[0]*unbindrad/(unbindrad-xtry[0])) : rate/(4*PI*difc*xtry[0]);
					errtry=(ratetry/rate-1.0)*(ratetry/rate-1.0)+(chitry/chi-1.0)*(chitry/chi-1.0);
					if(errtry<errold)												// better so go faster in same direction
						xmult[indx]*=1.1;
					else																		// worse so slow down and reverse
						xmult[indx]*=-0.5;
					if(errtry<errbest) {
						errbest=errtry;
						abest=xtry[0];
						pbest=xtry[1]; }}}}}

	if(chi>=0 && probptr)
		*probptr=pbest;

	return abest; }


/* unbindingradius returns the unbinding radius that corresponds to the geminate
reaction probability in pgem, the time step in dt, the mutual diffusion constant
in difc, and the binding radius in a.  Illegal inputs result in a return value
of -2.  If the geminate binding probability can be made as high as that
requested, the corresponding unbinding radius is returned.  Otherwise, the
negative of the maximum achievable pgem value is returned.
Modified 2/25/08 to allow for dt==0.  */
double unbindingradius(double pgem,double dt,double difc,double a) {
	double b,lo,dif,step,ki,kmax;
	int n;

	if(pgem<=0 || pgem>=1 || difc<=0 || a<=0 || dt<0) return -2;
	if(dt==0) return a/pgem;
	step=sqrt(2.0*difc*dt);
	ki=numrxnrate(step,a,-1);
	kmax=numrxnrate(step,a,0);
	if(1.0-ki/kmax<pgem) return ki/kmax-1.0;
	lo=0;
	b=a;
	while(1.0-ki/numrxnrate(step,a,b)>pgem) {
		lo=b;
		b*=2.0; }
	dif=b-lo;
	for(n=0;n<15;n++) {
		dif*=0.5;
		b=lo+dif;
		if(1.0-ki/numrxnrate(step,a,b)>pgem) lo=b; }
	b=lo+0.5*dif;
	return b; }


/******************************************************************************/
/*************    UTILITY FUNCTION, COPIED FROM MY MATH2.C FILE    ************/
/******************************************************************************/


double rxnparam_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; }


/******************************************************************************/
/************    FUNCTIONS FOR INVESTIGATING AN ABSORBING SPHERE    ***********/
/******************************************************************************/


/* rdfmodelrdf. Computes analytical model rdfs. */
double rdfmodelrdf(double r,double sigmab,double sigmau,double lambdap,double gamma) {
	double g,rp,sigmaup;

	if(lambdap<0 && gamma < 0) {				// Smoluchowski model
		if(sigmau<0) {
			if(r<=sigmab) g=0;
			else g=1-sigmab/r; }
		else {
			if(r<=sigmab) g=0;
			else if(r<sigmau) g=1-sigmab/r*(sigmau-r)/(sigmau-sigmab);
			else g=1; }}
	else if(gamma>=0) {								// Collins and Kimball model
		if(sigmau<0) {
			if(r<sigmab) g=0;
			else g=1-sigmab/r*sigmab/(sigmab+gamma); }
		else {
			if(r<sigmab) g=0;
			else if(r<sigmau) g=sigmab/r*sigmab*(sigmau-r)/(sigmau*gamma+sigmab*sigmau-sigmab*sigmab);
			else g=1; }}
	else {														// Doi model
			rp=r/sigmab;
			if(sigmau<0) {
				if(rp<=1) g=sinh(rp*lambdap)/(rp*lambdap*cosh(lambdap));
				else g=1-1/rp+tanh(lambdap)/(rp*lambdap); }
			else {
				sigmaup=sigmau/sigmab;
				if(rp<=1) g=sigmaup/rp*sinh(rp*lambdap)/sinh(lambdap)/(1+lambdap*(sigmaup-1)/tanh(lambdap));
				else if(rp<sigmaup) g=sigmaup/rp*(lambdap*(rp-1)*cosh(lambdap)+sinh(lambdap))/(lambdap*(sigmaup-1)*cosh(lambdap)+sinh(lambdap));
				else g=1; }}

	return g; }


/* rdfabsorb integrates the radial diffusion function (rdf) for 0<=r<=1, sets
those values to 0, and returns the integral.  r is a vector of radii; r[0] may
equal zero but that is not required; if not, then it is assumed that the rdf has
zero slope at the origin.
   Integration uses a spherical version of the trapezoid rule: at positions r0
and r1, the function f has values f0 and f1, leading to the linear interpolation
f=[(r-r0)f1+(r1-r)f0]/(r1-r0).  Its integral is A=Integral[4*pi*r^2*f(r)*dr]
A=4*pi/(r1-r0)*[(f1-f0)/4*(r1^4-r0^4)+(r1f0-r0f1)/3*(r1^3-r0^3)]
A=pi*(f1-f0)(r1+r0)(r1^2+r0^2)+4*pi/3*(r1f0-r0f1)(r1^2+r1r0+r0^2)
The left end of the integral assumes zero slope for the rdf.  The right end does
not terminate exactly at 1, but includes the upper left triangle of the final
trapezoid.  That way, if there are two absorptions in a row, the second one will
return an integral of 0, and area is properly conserved.  The problem is that it
does not terminate exactly at 1.  Furthermore, the correct relative location of
1 between two r[j] points depends on the function.  The best solution is to use
an unevenly spaced r[j] vector, with a very narrow separation about 1 and no
r[j] equal to 1.   */
double rdfabsorb(double *r,double *rdf,int n,double bindrad,double prob) {
	int j;
	double sum,r0,r1,f0,f1;

	r0=r1=0;
	f1=rdf[0];
	sum=0;
	for(j=(r[0]==0?1:0);r1<bindrad && j<n;j++) {
		r0=r1;
		f0=f1;
		r1=r[j];
		f1=rdf[j];
		sum+=PI*(f1-f0)*(r1+r0)*(r1*r1+r0*r0)+4.0*PI/3.0*(r1*f0-r0*f1)*(r1*r1+r1*r0+r0*r0); }
	f0=0;
	sum-=PI*(f1-f0)*(r1+r0)*(r1*r1+r0*r0)+4.0*PI/3.0*(r1*f0-r0*f1)*(r1*r1+r1*r0+r0*r0);
	sum*=prob;
	for(j-=2;j>=0;j--) rdf[j]*=(1.0-prob);
	return sum; }


/* rdfdiffuse integrates the radial distribution function with the Green's
function for radially symmetric diffusion to implement diffusion over a fixed
time step.  r is a vector of radii, rdfa is the input rdf, rdfd is the output
rdf, n is the number of points, and step is the rms step length, equal to
(2Dt)^1/2.  r[0] may equal 0 but it is not required.  It is assumed that rdfa
has zero slope at r=0.  The boundary condition on the large r side is that the
function tends to 1 with a functional form 1+a2/r, for large r.  This is
accomplished by fitting the 10% largest r portion of the rdf with the function
1+a2/r.  After the integral over the tabulated data is complete, the rdf is
integrated on to infinity using the previous fit information and an analytical
result for that integral.  The numerical portion of the integral is carried out
exactly like the one in absorb but with a different integrand, which is
c(r)=Integrate[4*pi*r'^2*rdfa(r')*grn(r,r')dr'].  grn(r,r') is the Green's function, equal
to grn(r,r')=1/(4*pi*r*r')[G_step(r-r')-G_step(r+r')] and G_s(x) is a normalized
Gaussian with mean 0 and standard deviation s. */
void rdfdiffuse(double *r,double *rdfa,double *rdfd,int n,double step) {
	int i,j;
	double grn,sum,f0,f1,rr,r0,r1,alpha,beta,a2,erfcdif,erfcsum;

	alpha=beta=0;												// fit the final 10% to 1+a2/r, finding a2 value
	for(i=(int)(0.9*n);i<n;i++) {
		alpha+=1.0/r[i]/r[i];
		beta+=(rdfa[i]-1.0)/r[i]; }
	a2=beta/alpha/step;

	grn=0;
	if(r[i=0]==0) {
		rr=r1=f1=sum=0;
		for(j=1;j<n;j++) {								// loop over r' values, for r=0
			r0=r1;
			f0=f1;
			r1=r[j]/step;
			grn=exp(-r1*r1/2.0)/(2.0*PI*SQRT2PI);		// green(0,r1)
			f1=(rdfa[j]-rdfa[0])*grn;
			sum+=PI*(f1-f0)*(r1+r0)*(r1*r1+r0*r0)+4.0*PI/3.0*(r1*f0-r0*f1)*(r1*r1+r1*r0+r0*r0); }
		sum+=(1.0-rdfa[0])*(4.0*PI*r1*grn+rxnparam_erfccD(r1/sqrt(2.0)));
		rdfd[i++]=sum+rdfa[0]; }
	for(;i<n;i++) {											// loop over r values
		rr=r[i]/step;
		r1=0;
		grn=exp(-rr*rr/2.0)/(2.0*PI*SQRT2PI);	// green(rr,0)
		f1=(rdfa[0]-rdfa[i])*grn;
		sum=0;
		for(j=(r[0]==0?1:0);j<n;j++) {		// loop over r' values
			r0=r1;
			f0=f1;
			r1=r[j]/step;
			grn=1.0/rr/r1*(exp(-(rr-r1)*(rr-r1)/2.0)-exp(-(rr+r1)*(rr+r1)/2.0))/(4.0*PI*SQRT2PI);
			f1=(rdfa[j]-rdfa[i])*grn;
			sum+=PI*(f1-f0)*(r1+r0)*(r1*r1+r0*r0)+4.0*PI/3.0*(r1*f0-r0*f1)*(r1*r1+r1*r0+r0*r0); }
		erfcdif=rxnparam_erfccD((r1-rr)/sqrt(2.0));
		erfcsum=rxnparam_erfccD((r1+rr)/sqrt(2.0));
		sum+=(1.0-rdfa[i])*(4.0*PI*r1*grn+1.0/2.0*(erfcdif+erfcsum))+a2/2.0/rr*(erfcdif-erfcsum);
		rdfd[i]=sum+rdfa[i]; }
	return; }


/* Analysis of the reverse reaction involves adding a delta function to the rdf
and then convolving with the Green's function.  However, this leads to
numerical errors, although it is trivial analytically.  The function
rdfreverserxn is the analytic solution.  It adds the diffusion Green's function
to the rdf, based on a delta function at b, and after one diffusion step.  r is
a list of radii, rdf is the rdf, step is the rms step length, b is the delta
function point (which does not have to be equal or unequal to a r[j] value), and
flux is the area of the delta function. */
void rdfreverserxn(double *r,double *rdf,int n,double step,double b,double flux) {
	int i;
	double rr,k;

	k=1.0/(4.0*PI*SQRT2PI*step*step*step);		// the 4*PI*SQRT2PI are from the Green's function
	if(b==0) {
		for(i=0;i<n;i++) {
			rr=r[i]/step;
			rdf[i]+=flux*k*2.0*exp(-rr*rr/2.0); }}
	else {
		b/=step;
		if(r[i=0]==0) rdf[i++]+=flux*k*2.0*exp(-b*b/2.0);
		for(;i<n;i++) {
			rr=r[i]/step;
			rdf[i]+=flux*k/rr/b*(exp(-(rr-b)*(rr-b)/2.0)-exp(-(rr+b)*(rr+b)/2.0)); }}
	return; }


/* rdfsteadystate calculates the radial distribution function (rdf) for alternating
absorption and diffusion steps, for either irreversible or reversible reactions.
r is a vector of radii, rdfa is input as a trial rdf and output as the result
after absorption, rdfd is ignored on input but is output as the rdf after
diffusion, n is the number of elements in the vectors, step is the rms step
length, and b is either <0 if the reaction is irreversible or is the unbinding
radius if the reaction is reversible.  It executes until the fractional
difference between successive steps is less than eps, but at least 30 times and
no more than maxit times.  It can also exit if it iterates more than maxit times
before converging or if the flux exceeds maxflux; if either of these happens,
the function returns -1. */
double rdfsteadystate(double *r,double *rdfa,double *rdfd,int n,double step,double a,double b,double eps,double prob) {
	const int maxit=100000;
	const double maxflux=1e7;
	int i,it;
	double flux,fluxold;

	rdfdiffuse(r,rdfa,rdfd,n,step);
	flux=fluxold=rdfabsorb(r,rdfd,n,a,prob);
	for(it=0;it<30 || (it<maxit && flux<maxflux && fabs((flux-fluxold)/(fluxold+1e-20))>eps);it++) {
		fluxold=flux;
		rdfdiffuse(r,rdfa,rdfd,n,step);
		if(b>=0) rdfreverserxn(r,rdfd,n,step,b,flux);
		for(i=0;i<n;i++) rdfa[i]=rdfd[i];
		flux=rdfabsorb(r,rdfa,n,a,prob); }
	if(it>=maxit || flux>=maxflux) flux=-1;
	return flux; }


/* rdfmaketable is used to create data tables of reaction rates, including those
used above in numrxnrate.  All input is requested from the user using the
standard input and all output is sent to standard output.
Runtime is about 1 minute with mode i, 200 pts, eps=1e-4 */
void rdfmaketable() {
	double slo=exp(-3.0),shi=exp(3.0),sinc=exp(0.2);		// step size low, high, increment
	const double blor=exp(0.0),bhir=exp(3.0),bincr=exp(0.2);	// b value low, high, increment for b>a
	const double blob=0,bhib=1.0,bincb=0.1;					// b value low, high, increment for b<a
	double *r,*rdfa,*rdfd,dr,s,b,flux,eps;
	int i,n,done;
	char mode,dir,string[256];

	printf("Function for calculating radial diffusion functions (rdf) and reactive\n");
	printf("fluxes for alternating reaction and diffusion steps.  This module\n");
	printf("can operate in several modes.  Enter (i) for irreversible reactions\n");
	printf("(r) for reversible reactions with the unbinding radius larger than\n");
	printf("the binding radius, or (b) for other reversible reactions.  Enter this\n");
	printf("mode in upper case for machine readable output.\n");
	printf("Operation mode: ");
	scanf("%s",string);
	mode=string[0];
	printf("Enter the number of radial points in the rdf (e.g. 200): ");
	scanf("%i",&n);
	if(n<10) {
		printf("Value is too low.  Function stopped.\n");return; }
	printf("Enter level of precision (e.g. 1e-4): ");
	scanf("%lf",&eps);
	if(eps<=0) {
		printf("Impossible precision.  Function stopped.\n");return; }
	printf("Enter u for increasing step lengths, d for decreasing: ");
	scanf("%s",string);
	dir=string[0];
	if(dir=='d') {
		s=slo;slo=shi;shi=s;
		sinc=1.0/sinc; }

	r=(double*)calloc(n,sizeof(double));
	rdfa=(double*)calloc(n,sizeof(double));
	rdfd=(double*)calloc(n,sizeof(double));
	if(!r || !rdfa || !rdfd) {
		printf("Out of memory.  Function stopped.\n");return; }

	if(mode=='i' || mode=='I') b=-1;
	else if(mode=='r' || mode=='R') b=blor;
	else b=blob;
	done=0;
	if(mode=='i') printf("step     flux\n");
	else if(mode=='r' || mode=='b') printf("b      step       flux\n");
	else printf("\n");

	while(!done) {
		if(mode=='i' || mode=='I') dr=10.0/n;						// r max is set to 10 for mode i
		else if(mode=='r' || mode=='R') dr=(b+3.0)/n;		// r max is set to b+3 for mode r
		else dr=5.0/n;																// r max is set to 5 for mode b
		r[0]=0;																				// r min is set to 0
		for(i=1;r[i-1]<1 && i-1<n;i++) r[i]=r[i-1]+dr;
		r[i-1]=0.9999;																// r point below 1 is set to 0.9999
		r[i++]=1.0001;																// r point above 1 is set to 1.0001
		for(;i<n;i++) r[i]=r[i-1]+dr;

		for(i=0;i<n && r[i]<1;i++) rdfa[i]=0;
		if(dir=='u' && (mode=='i' || mode=='I'))
			for(;i<n;i++) rdfa[i]=1.0-1.0/r[i];
		else if(dir=='u' && (mode=='r' || mode=='R'))
			for(;i<n && r[i]<b;i++) rdfa[i]=1.0-(b-r[i])/r[i]/(b-1.0);
		for(;i<n;i++) rdfa[i]=1.0;

		flux=0;
		for(s=slo;slo<shi?s<shi*sqrt(sinc):s>shi*sqrt(sinc);s*=sinc) {
			if(flux>=0) flux=rdfsteadystate(r,rdfa,rdfd,n,s,1,b,eps,1);
			if(mode=='i') printf("%lf %lf\n",s,flux);
			else if(mode=='r' || mode=='b') printf("%lf %lf %lf\n",b,s,flux);
			else printf("%lf,",fabs(flux)); }
		printf("\n");

		if(mode=='i' || mode=='I') done=1;
		else if(mode=='r' || mode=='R') {
			b*=bincr;
			if(b>bhir*sqrt(bincr)) done=1; }
		else {
			b+=bincb;
			if(b>bhib+bincb/2.0) done=1; }}
	free(r);
	free(rdfa);
	free(rdfd);
	return; }


/* rdfmaketableprob is identical to rdfmaketable but also includes a reaction
probability. It also uses slightly different default parameters. */
void rdfmaketableprob() {
	double slo=exp(-3.0),shi=exp(3.0),sinc=exp(0.2);		// step size low, high, increment
	const double blor=exp(0.0),bhir=exp(3.0),bincr=exp(0.2);	// for b>=1; b value low, high, increment for b>a
	const double blob=0,bhib=1.0,bincb=0.1;					// for b<=1; b value low, high, increment for b<a
	const double pvect[]={0,0.01,0.02,0.05,0.1,0.2,0.5,1.0};
	const int pnum=8;
	double *r,*rdfa,*rdfd,dr,s,b,prob,flux,eps;
	int i,n,done,it,pindx;
	char mode,dir,string[256];

	printf("\nFunction for calculating radial diffusion functions (rdf) and reactive\n");
	printf("fluxes for alternating reaction and diffusion steps.  This module\n");
	printf("can operate in several modes.  Enter (i) for irreversible reactions\n");
	printf("(r) for reversible reactions with the unbinding radius larger than\n");
	printf("the binding radius, or (b) for reversible reactions with unbinding\n");
	printf("inside the binding radius.  Enter this mode in upper case for machine\n");
	printf("readable output.\n");
	printf("Operation mode: ");
	scanf("%s",string);
	mode=string[0];
	printf("Enter the number of radial points in the rdf (e.g. 200): ");
	scanf("%i",&n);
	if(n<10) {
		printf("Value is too low.  Function stopped.\n");return; }
	printf("Enter level of precision (e.g. 1e-6): ");
	scanf("%lf",&eps);
	if(eps<=0) {
		printf("Impossible precision.  Function stopped.\n");return; }
	if(mode=='b'||mode=='B') {
		printf("Using decreasing step lengths because in mode b or B.\n");
		dir='d'; }
	else {
		printf("Enter u for increasing step lengths, d for decreasing: ");
		scanf("%s",string);
		dir=string[0]; }
	if(dir=='d') {
		s=slo;slo=shi;shi=s;
		sinc=1.0/sinc; }
	printf("\n");

	r=(double*)calloc(n,sizeof(double));
	rdfa=(double*)calloc(n,sizeof(double));
	rdfd=(double*)calloc(n,sizeof(double));
	if(!r || !rdfa || !rdfd) {
		printf("Out of memory.  Function stopped.\n");return; }

	if(mode=='i' || mode=='I') b=-1;									// i/I mode is for irreversible.  I for machine readable.
	else if(mode=='r' || mode=='R') b=blor;						// r,R mode is for reversible with b>=1.  R for machine readable.
	else b=blob;																			// b,B mode is for reversible with b<=1.  B for machine readable.
	done=0;
	if(mode=='i') printf("iter  prob    step   flux\n");
	else if(mode=='r' || mode=='b') printf("iter  sigma_u    prob     step       flux\n");
	else {
		if(mode=='R')
			printf("b from %g to %g with steps %g\n",blor,bhir,bincr);
		else if(mode=='B')
			printf("b from %g to %g with steps %g\n",blob,bhib,bincb);
		printf("prob from %g to %g in %i steps\n",pvect[0],pvect[pnum-1],pnum);
		printf("step from %g to %g with steps %g\n\n",slo,shi,sinc); }

	it=0;
	while(!done) {
		if(mode=='i' || mode=='I') dr=10.0/n;						// r max is set to 10 for mode i
		else if(mode=='r' || mode=='R') dr=(b+3.0)/n;		// r max is set to b+3 for mode r
		else dr=5.0/n;																	// r max is set to 5 for mode b
		r[0]=0;																					// r min is set to 0
		for(i=1;r[i-1]<1 && i-1<n;i++) r[i]=r[i-1]+dr;
		r[i-1]=0.9999;																	// r point below 1 is set to 0.9999
		r[i++]=1.0001;																	// r point above 1 is set to 1.0001
		for(;i<n;i++) r[i]=r[i-1]+dr;

		for(pindx=0;pindx<pnum;pindx++) {
			prob=pvect[pindx];

			for(i=0;i<n && r[i]<1;i++) rdfa[i]=0;
			if(dir=='u' && (mode=='i' || mode=='I'))
				for(;i<n;i++) rdfa[i]=1.0-1.0/r[i];					// use of Smoluchowski model is not ideal but I'm not sure what's better
			else if(dir=='u' && (mode=='r' || mode=='R'))
				for(;i<n && r[i]<b;i++) rdfa[i]=1.0-(b-r[i])/r[i]/(b-1.0);
			for(;i<n;i++) rdfa[i]=1.0;

			flux=0;
			for(s=slo;slo<shi?s<shi*sqrt(sinc):s>shi*sqrt(sinc);s*=sinc) {
				it++;
				if(prob==0) flux=0;
				else if(flux>=0) flux=rdfsteadystate(r,rdfa,rdfd,n,s,1,b,eps,prob);
				if(mode=='i') printf("%i %lf %lf %lf\n",it,prob,s,flux);
				else if(mode=='r' || mode=='b') printf("%i %lf %lf %lf %lf\n",it,b,prob,s,flux);
				else printf("%lf,",flux); }
			printf("\n"); }

		if(mode=='i' || mode=='I') done=1;
		else if(mode=='r' || mode=='R') {
			b*=bincr;
			if(b>bhir*sqrt(bincr)) done=1; }
		else {
			b+=bincb;
			if(b>bhib+bincb/2.0) done=1; }}
	free(r);
	free(rdfa);
	free(rdfd);
	return; }