/* Steven Andrews, 11/20/03 */
/* Rewrite of rebind.c, because that program is confusing. */
/* Here, either of two systems are simulated: clustered receptors on a sphere
and unclustered receptors on a sphere.  In the former case, the receptors extend
over only a small portion of the sphere; in the latter case, they cover the
whole sphere nearly uniformly. */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "random.h"

#define PI 3.14159265358979323846

#define DIFF 500.0						// ### diffusion constant, in nm^2/us  (500 for aspartate in water, 136 CheR in cyto)
#define TAUU 1000							// ### unbinding time constant, in us (1000 for Tar-aspartate, 45000 for Tar-CheR)
#define AVALUE 0.476					// ### a value, in nm (1.0 default, 0.264 aspartate, 0.0019 CheR)
#define BVALUE 4.756					// ### b value, in nm (3.0 default, 0.528 aspartate, 0.19 CheR)
#define DXMIN (AVALUE/100.0)	// ### minimum rms step length, in nm
#define SPHRAD 750.0					// ### radius of sphere, in nm (750 for E coli)
#define NRECEPT 3000					// ### number of receptors (3000 for E coli)
#define CLUSTERED 1						// &&& flag for clustered receptors (1) vs. unclustered (0)
#define CLRAD 225.0						// ### radius of cluster by arc length (225 default); only used if clustered
#define ESCAPE (SPHRAD*1000.0)	// ### escape radius
#define RANDTABLEMAX 4095			// ### size of random number table minus 1

float GAUSS[RANDTABLEMAX+1];	// table of random numbers

void diffuse(double *x,double dx);
double spherelatticedist(double *xptr,double *yptr,double *zptr,int n,double rad);
double spherelatticedist2(double *xptr,double *yptr,double *zptr,int n,double rad,double latmin);
void spherelatticepts(int n,double rad,double latmin,FILE *fptr);
void spherelatticerows(int n,double rad,double latmin,FILE *fptr);

double dist2recept(double x[3],double xr[3]);
double dist2surf(double x[3]);
void reflectsurface(double x[3]);
void move2unbind(double x[3]);
int diffuse2bind(double x[3],double *tptr);

void geminatebind(void);
void initialsite(void);
void allsites(void);
void radcorrfn(void);
void radcorrall(void);
void barfigure(void);
void bindrate(void);
void equilibrium(void);
void outputrecept(void);
void teststuff(void);


/* ********************* GENERAL PURPOSE ROUTINES *********************** */


/* diffuse carries out 3-D diffusion on vector x, using rms step length dx. */
void diffuse(double *x,double dx) {
	x[0]+=dx*GAUSS[rand()&RANDTABLEMAX];
	x[1]+=dx*GAUSS[rand()&RANDTABLEMAX];
	x[2]+=dx*GAUSS[rand()&RANDTABLEMAX];
	return; }


/* spherelatticedist returns the distance to the nearest lattice site that is on
the surface of a sphere.  The sphere is centered about the origin and has radius
rad.  The lattice on the surface has about n points, although there is some
variation due to discretization effects.  Send in the current position in xptr,
yptr, and zptr.  These are returned with the location of the nearest lattice
site.  The value returned by the function is the distance to that point.
Lattice sites are arranged in rows along equally spaced latitudes.  The lattice
is square near the equator and prime meridian for n large, and is non-square
elsewhere. */
double spherelatticedist(double *xptr,double *yptr,double *zptr,int n,double rad) {
	double x,y,z,r,lat,phi;
	double deltalat,lat2,deltaphi,phi2;
	int latn,phin,lati;

	x=*xptr;
	y=*yptr;
	z=*zptr;
	r=sqrt(x*x+y*y+z*z);
	lat=(r==0)?0:asin(z/r);
	phi=atan2(y,x);
	latn=floor(sqrt(PI*n)/2.0+0.5);										// number of rows
	deltalat=PI/latn;																	// spacing of rows
	if(latn%2) {
		lati=floor(lat/deltalat+0.5);										// row number
		if(lati>latn/2) lati=latn/2;
		else if(lati<-latn/2) lati=-latn/2;
		lat2=lati*deltalat; }
	else {
		lati=floor(lat/deltalat);
		if(lati>=latn/2) lati=latn/2-1;
		else if(lati<-latn/2) lati=-latn/2;
		lat2=(lati+0.5)*deltalat; }
	phin=floor(sqrt(n*PI)*cos(lat2)+0.5);							// number of points on row
	deltaphi=2*PI/phin;																// spacing along row
	phi2=floor(phi/deltaphi+0.5)*deltaphi;						// long of point
	*xptr=rad*cos(lat2)*cos(phi2);
	*yptr=rad*cos(lat2)*sin(phi2);
	*zptr=rad*sin(lat2);
	return sqrt((x-*xptr)*(x-*xptr)+(y-*yptr)*(y-*yptr)+(z-*zptr)*(z-*zptr)); }


/* spherelatticedist2 returns the distance to the nearest lattice site that is
on the surface of a sphere and that is "north" of the latitute latmin.  The
sphere is centered about the origin and has radius rad.  The lattice within the
patch above latmin has about n points, although there is some variation due to
discretization effects.  Send in the current position in xptr, yptr, and zptr.
These are returned with the location of the nearest lattice site.  The value
returned by the function is the distance to that point in a straight line
distance.  Lattice sites are arranged in rows along equally spaced latitudes.
The lattice is square near the equator and prime meridian for n large, and is
non-square elsewhere. */
double spherelatticedist2(double *xptr,double *yptr,double *zptr,int n,double rad,double latmin) {
	double x,y,z,r,lat,phi,f;
	double deltalat,lat2,deltaphi,phi2;
	int latn,phin,lati;

	x=*xptr;
	y=*yptr;
	z=*zptr;
	r=sqrt(x*x+y*y+z*z);
	lat=(r==0)?0:asin(z/r);
	phi=atan2(y,x);
	f=(1.0-sin(latmin))/2.0;
	latn=floor(sqrt(PI*n/f)/2.0+0.5);									// number of rows
	deltalat=PI/latn;																	// spacing of rows
	if(lat<latmin) lat=latmin;
	if(latn%2) {
		lati=floor(lat/deltalat+0.5);										// row number
		if(lati>latn/2) lati=latn/2;
		else if(lati<-latn/2) lati=-latn/2;
		lat2=lati*deltalat; }
	else {
		lati=floor(lat/deltalat);
		if(lati>=latn/2) lati=latn/2-1;
		else if(lati<-latn/2) lati=-latn/2;
		lat2=(lati+0.5)*deltalat; }
	phin=floor(sqrt(n/f*PI)*cos(lat2)+0.5);						// number of points on row
	deltaphi=2*PI/phin;																// spacing along row
	phi2=floor(phi/deltaphi+0.5)*deltaphi;						// long of point
	*xptr=rad*cos(lat2)*cos(phi2);
	*yptr=rad*cos(lat2)*sin(phi2);
	*zptr=rad*sin(lat2);
	return sqrt((x-*xptr)*(x-*xptr)+(y-*yptr)*(y-*yptr)+(z-*zptr)*(z-*zptr)); }


/* spherelattticepts lists all points for a lattice on a sphere.  As usual, n is
the approximate number of points, rad is the sphere radius, and latmin is the
minimum latitude on the sphere where there are points.  fptr is the file stream
for the output.  These points exactly correspond to those of the other routines. */
void spherelatticepts(int n,double rad,double latmin,FILE *fptr) {
	int latn,lati,phin,phii,npts;
	double f,deltalat,lat2,deltaphi,phi2,x,y,z;

	npts=0;
	f=(1.0-sin(latmin))/2.0;
	latn=floor(sqrt(PI*n/f)/2.0+0.5);									// number of rows
	deltalat=PI/latn;																	// spacing of rows
	if(latn%2) {
		lati=floor(latmin/deltalat+0.5);
		if(lati>latn/2) lati=latn/2;
		else if(lati<-latn/2) lati=-latn/2; }
	else {
		lati=floor(latmin/deltalat);
		if(lati>=latn/2) lati=latn/2-1;
		else if(lati<-latn/2) lati=-latn/2; }
	for(;lati<=((latn%2)?latn/2:latn/2-1);lati++) {
		lat2=(latn%2)?lati*deltalat:(lati+0.5)*deltalat;
		phin=floor(sqrt(n/f*PI)*cos(lat2)+0.5);						// number of points on row
		deltaphi=2*PI/phin;																// spacing along row
		for(phii=0;phii<phin;phii++) {
			phi2=phii*deltaphi;
			x=rad*cos(lat2)*cos(phi2);
			y=rad*cos(lat2)*sin(phi2);
			z=rad*sin(lat2);
			if(fptr) fprintf(fptr,"%lf %lf %lf\n",x,y,z);
			else printf("%lf %lf %lf\n",x,y,z);
			npts++; }}
	printf("n=%i\n",npts);
	return; }


/* spherelattticerows lists all rows for a lattice on a sphere.  As usual, n is
the approximate number of points, rad is the sphere radius, and latmin is the
minimum latitude on the sphere where there are points.  fptr is the file stream
for the output.  These points exactly correspond to those of the other routines. */
void spherelatticerows(int n,double rad,double latmin,FILE *fptr) {
	int latn,lati,phin,phii,npts;
	double f,deltalat,lat2,deltaphi,phi2,x,y,z;

	npts=0;
	f=(1.0-sin(latmin))/2.0;
	latn=floor(sqrt(PI*n/f)/2.0+0.5);									// number of rows
	deltalat=PI/latn;																	// spacing of rows
	if(latn%2) {
		lati=floor(latmin/deltalat+0.5);
		if(lati>latn/2) lati=latn/2;
		else if(lati<-latn/2) lati=-latn/2; }
	else {
		lati=floor(latmin/deltalat);
		if(lati>=latn/2) lati=latn/2-1;
		else if(lati<-latn/2) lati=-latn/2; }
	printf("lat  phin  deltaphi  phisp\n");
	for(;lati<=((latn%2)?latn/2:latn/2-1);lati++) {
		lat2=(latn%2)?lati*deltalat:(lati+0.5)*deltalat;
		phin=floor(sqrt(n/f*PI)*cos(lat2)+0.5);						// number of points on row
		deltaphi=2*PI/phin;																// spacing along row
		if(fptr) fprintf(fptr,"%lf %i\n",lat2,phin);
		else printf("%lf %i %lf %lf\n",lat2,phin,deltaphi,deltaphi*rad*cos(lat2));
		npts+=phin; }
	printf("latmin=%lf f=%lf\n",latmin,f);
	printf("latn=%i deltalat=%lf latsp=%lf\n",latn,deltalat,rad*deltalat);
	printf("n: %i\n",npts);
	return; }







/* ****************** ROUTINES SPECIFIC TO THIS PROGRAM ***************** */


/* Returns distance to nearest receptor and, in xr, the position of that
receptor.  To speed this routine up, enter the position of the previous closest
receptor in xr and if it is still probably the closest one, then the receptor
isn't identified again. */
double dist2recept(double x[3],double xr[3]) {
	double r,f,deltalat;
	int latn;
	static double rmin=0;

	if(rmin==0) {
		if(CLUSTERED) f=(1.0-sin(PI/2.0-CLRAD/SPHRAD))/2.0;
		else f=1.0;
		latn=floor(sqrt(PI*NRECEPT/f)/2.0+0.5);						// number of rows
		deltalat=PI/latn;																	// spacing of rows
		rmin=SPHRAD*deltalat*2.0/3.0; }
	r=sqrt((x[0]-xr[0])*(x[0]-xr[0])+(x[1]-xr[1])*(x[1]-xr[1])+(x[2]-xr[2])*(x[2]-xr[2]));
	if(r<rmin) return r;
	xr[0]=x[0];
	xr[1]=x[1];
	xr[2]=x[2];
	#if CLUSTERED==0
		r=spherelatticedist(&xr[0],&xr[1],&xr[2],NRECEPT,SPHRAD);
	#else
		r=spherelatticedist2(&xr[0],&xr[1],&xr[2],NRECEPT,SPHRAD,PI/2.0-CLRAD/SPHRAD);
	#endif
	return r; }


/* dist2surf returns the distance from x to the surface of the sphere.  If x is
inside the sphere, it returns a negative number. */
double dist2surf(double x[3]) {
	return sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])-SPHRAD; }


/* Executes reflection off the surface of a sphere, assuming the distance to the
sphere surface is small, so that errors are negligible. */
void reflectsurface(double x[3]) {
	double d;

	d=SPHRAD*SPHRAD/(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
	x[0]*=d;
	x[1]*=d;
	x[2]*=d;
	return; }


/* move2unbind inputs a point that is a receptor point (this is not checked) and
replaces it with the unbinding location for that receptor. */
void move2unbind(double x[3]) {
	double d;

	d=(SPHRAD+BVALUE)/SPHRAD;
	x[0]*=d;
	x[1]*=d;
	x[2]*=d;
	return; }


/* diffuse2bind diffuses a ligand from x until it either binds a receptor or
escapes the system.  If it binds a receptor, the function returns 1, otherwise
it returns zero.  The returned x value is the receptor location if it ends up
bound.  tptr is a pointer to the time variable; it is returned with the time
incremented for the amount of time that the ligand diffused.  */
int diffuse2bind(double x[3],double *tptr) {
	int bind,escape;
	double dx,xr[3],d2m,d2s,d2r,rspc;

	rspc=sqrt(4.0*PI*SPHRAD*SPHRAD/NRECEPT);
	dx=DXMIN;
	xr[0]=xr[1]=xr[2]=0;
	bind=escape=0;
	while(!bind&&!escape) {
		diffuse(x,dx);
		*tptr+=dx*dx/(2.0*DIFF);
		d2m=sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
		d2s=d2m-SPHRAD;
		if(d2s<rspc) d2r=dist2recept(x,xr)-AVALUE;
		else d2r=d2s;
		if(d2r<=0) bind=1;
		else if(d2s<0) reflectsurface(x);
		else if(d2m>ESCAPE) escape=1;
		dx=0.1*(1.0-0.9*SPHRAD/d2m)*d2m;
		if(d2r<rspc&&0.2*d2r<dx) dx=0.2*d2r;
		if(dx<DXMIN) dx=DXMIN; }
	if(bind) {
		x[0]=xr[0];
		x[1]=xr[1];
		x[2]=xr[2]; }
	return bind; }






/* **************** TOP LEVEL ROUTINES ********************** */





/* finds geminate rebinding probability using model system */
void geminatebind(void) {
	double x[3],xi[3];					// location of ligand, initial receptor
	int itmax=100000;						// ### number of iterations
	int it,gem;									// iterations and number of geminate bindings
	double t;										// time (ignored)
	int bind;										// flag for binding

	gem=0;
	for(it=0;it<itmax;it++) {
		t=0;
		sphererandd(x,SPHRAD,SPHRAD);
		xi[0]=xi[1]=xi[2]=0;
		dist2recept(x,xi);
		x[0]=xi[0];
		x[1]=xi[1];
		x[2]=xi[2];
		move2unbind(x);						// starting position
		bind=diffuse2bind(x,&t);
		if(bind&&x[0]==xi[0]&&x[1]==xi[1]&&x[2]==xi[2]) gem++;
		if(!((it+1)%1000))
			printf("gem, total, ratio: %i %i %lf\n",gem,it+1,(double)gem/(it+1)); }
	return; }


/* equilibrium finds the ratio of time that a ligand spends bound to a receptor
to the time that it is not bound to a receptor.  The system here is the standard
model system except that it is enclosed in a smallish box with periodic boundary
conditions.  Simulation time is 190 seconds for tmax=100000000 */
void equilibrium(void) {
	double x[3],xr[3],dx,d2s,d2m,d2r,rspc;
	double xlo,xhi,ylo,yhi,zlo,zhi;
	double t,ton,tont;					// time, time on, time on total
	double tmax=100000000;			// ### maximum time
	int bind,it;
	double vol,ratio;

	it=0;
	xlo=ylo=zlo=-1.5*SPHRAD;		// ### outside of box
	xhi=yhi=zhi=1.5*SPHRAD;			// ### outside of box
	vol=(xhi-xlo)*(yhi-ylo)*(zhi-zlo)-4.0/3.0*PI*SPHRAD*SPHRAD*SPHRAD;
	printf("volume: %e\n",vol);
	ratio=2.0*PI*DIFF*AVALUE*NRECEPT*TAUU/((1-AVALUE/BVALUE)*vol);
	printf("correct on/off ratio: %f\n",ratio);
	sphererandd(x,1.1*SPHRAD,1.1*SPHRAD);
	xr[0]=xr[1]=xr[2]=0;
	bind=0;
	tont=0;
	rspc=sqrt(4.0*PI*SPHRAD*SPHRAD/NRECEPT);
	for(t=0;t<tmax;) {
		dx=DXMIN;
		while(!bind&&t<tmax) {
			diffuse(x,dx);
			t+=dx*dx/(2.0*DIFF);
			d2m=sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
			d2s=d2m-SPHRAD;
			if(d2s<rspc) d2r=dist2recept(x,xr)-AVALUE;
			else d2r=d2s;
			if(d2r<=0) bind=1;
			else if(d2s<0) reflectsurface(x);
			else if(x[0]<xlo) x[0]=xhi-(xlo-x[0]);
			else if(x[0]>xhi) x[0]=xlo+(x[0]-xhi);
			else if(x[1]<ylo) x[1]=yhi-(ylo-x[1]);
			else if(x[1]>yhi) x[1]=ylo+(x[1]-yhi);
			else if(x[2]<zlo) x[2]=zhi-(zlo-x[2]);
			else if(x[2]>zhi) x[2]=zlo+(x[2]-zhi);
			dx=0.1*(1.0-0.9*SPHRAD/d2m)*d2m;
			if(d2r<rspc&&0.2*d2r<dx) dx=0.2*d2r;
			if(dx<DXMIN) dx=DXMIN; }
		if(bind) {
			x[0]=xr[0];
			x[1]=xr[1];
			x[2]=xr[2];
			move2unbind(x);
			ton=exprand30(TAUU);
			tont+=ton;
			t+=ton;
			bind=0;
			it++; }
		if(!(it%100)||t>=tmax) printf("bindings, time, time on / time off: %i %f %f / %f = %f\n",it,t,tont,t-tont,tont/(t-tont)); }
	return; }


/* finds initial binding sites of ligands that start far away.  Runtime is 10
minutes for itmax=10000. */
void initialsite(void) {
	double x[3],t;							// location of ligand, time
	int itmax=10000;						// ### number of iterations
	int it,bind;								// iteration, bind flag
	FILE *fptr;

	fptr=fopen("initsite","w");
	t=0;
	for(it=0;it<itmax;) {
		sphererandd(x,1.1*SPHRAD,1.1*SPHRAD);
		bind=diffuse2bind(x,&t);
		if(bind) {
			fprintf(fptr,"%lf\n",SPHRAD*asin(sqrt(x[0]*x[0]+x[1]*x[1])/SPHRAD));
			it++; }}
		fclose(fptr);
	return; }


/* finds all binding sites of ligands that start far away.  Runtime is 75
seconds for itmax=10000. */
void allsites(void) {
	double x[3],xo[3],t;							// location of ligand, time
	int itmax=10000;						// ### number of iterations
	int it,bind,sum,gem;						// iteration, bind flag
	FILE *fptr;

	fptr=fopen("allsite","w");
	t=0;
	bind=0;
	sum=0;
	for(it=0;it<itmax;) {
		if(!bind) sphererandd(x,1.1*SPHRAD,1.1*SPHRAD);
		bind=diffuse2bind(x,&t);
		if(bind) {
			sum++;
			gem=sum>1&&x[0]==xo[0]&&x[1]==xo[1]&&x[2]==xo[2];
			fprintf(fptr,"%i %i %lf\n",sum==1,gem,SPHRAD*asin(sqrt(x[0]*x[0]+x[1]*x[1])/SPHRAD));
			xo[0]=x[0];xo[1]=x[1];xo[2]=x[2];
			move2unbind(x);
			it++; }
		else sum=0; }
		fclose(fptr);
	return; }


/* finds binding sites of ligand that start at the center.  Runtime is 43
seconds for itmax=10000. */
void radcorrfn(void) {
	double x[3],xo[3],t,th;			// location of ligand, time
	int itmax=10000;						// ### number of iterations
	int it,bind,esc;							// iteration, bind flag
	FILE *fptr;

	fptr=fopen("radcorrfn","w");
	t=0;
	esc=0;
	for(it=0;it<itmax;) {
		th=unirand(0,2*PI);
		x[0]=1.0*cos(th);
		x[1]=1.0*sin(th);
		x[2]=SPHRAD;
		spherelatticedist2(&x[0],&x[1],&x[2],NRECEPT,SPHRAD,PI/2.0-CLRAD/SPHRAD);
		xo[0]=x[0];
		xo[1]=x[1];
		xo[2]=x[2];
		move2unbind(x);
		bind=diffuse2bind(x,&t);
		if(bind&&(x[0]!=xo[0]||x[1]!=xo[1]||x[2]!=xo[2])) {
			fprintf(fptr,"%f\n",sqrt((x[0]-xo[0])*(x[0]-xo[0])+(x[1]-xo[1])*(x[1]-xo[1])));
			it++; }
		else esc++; }
		fprintf(fptr,"%i\n",esc);
		fclose(fptr);
	return; }


/* finds all binding sites of ligand that start at the center.  Runtime is 23
seconds for itmax=10000. */
void radcorrall(void) {
	double x[3],t;			// location of ligand, time
	int itmax=10000;						// ### number of iterations
	int it,bind,esc;							// iteration, bind flag
	FILE *fptr;

	fptr=fopen("radcorrfn","w");
	t=0;
	bind=esc=0;
	for(it=0;it<itmax;) {
		if(!bind) {
			x[0]=x[1]=0;
			x[2]=SPHRAD; }
		bind=diffuse2bind(x,&t);
		if(bind) {
			fprintf(fptr,"%f\n",SPHRAD*asin(sqrt(x[0]*x[0]+x[1]*x[1])/SPHRAD));
			move2unbind(x);
			it++; }
		else esc++; }
		fprintf(fptr,"%i\n",esc);
		fclose(fptr);
	return; }


/* barfigure creates a postscript file with a stripe for each trial.  Black
sections are for times when a ligand is bound and white are for when it's not
bound. */

void barfigure(void) {
	double x[3],t,to,sf;
	double tmax=100;
	double mingap=0.00001;
	int row,bind;
	int rowmin=1;
	int rowmax=10;
	int seed=0;
	FILE *fptr;

	seed+=CLUSTERED*100;
	fptr=fopen("barfigure","w");
	fprintf(fptr,"%%!\n");
	fprintf(fptr,"%%%% Postscript file\n");
	sf=7.0/tmax;											// first number is final width in inches
	fprintf(fptr,"%%%% scaling factor: %f\n",sf);
	fprintf(fptr,"/inch {72 mul} def\n");
	fprintf(fptr,"3 setlinewidth\n");
	fprintf(fptr,"0.7 setgray\n");
	for(row=rowmin;row<=rowmax;row++) {
		fprintf(fptr,"newpath\n");
		fprintf(fptr,"%f inch %f inch moveto\n",0.5,5.0-0.1*row);
		fprintf(fptr,"%f inch %f inch lineto\n",0.5+tmax*sf,5.0-0.1*row);
		fprintf(fptr,"stroke\n"); }
	fprintf(fptr,"0 setgray\n");
	for(row=rowmin;row<=rowmax;) {
		srand(seed++);
		sphererandd(x,1.1*SPHRAD,1.1*SPHRAD);
		bind=diffuse2bind(x,&t);
		if(bind) {
			fprintf(fptr,"newpath\n");
			for(t=to=0;t<tmax&&bind;) {
				if(t!=to&&(t-to)*sf<mingap) t=to+mingap/sf;
				fprintf(fptr,"%f inch %f inch moveto\n",0.5+t*sf,5.0-0.1*row);
				t+=exprand(TAUU);
				move2unbind(x);
				if(t>tmax) t=tmax;
				fprintf(fptr,"%f inch %f inch lineto\n",0.5+t*sf,5.0-0.1*row);
				to=t;
				bind=diffuse2bind(x,&t); }
			fprintf(fptr,"stroke\n");
			row++; }}
	fprintf(fptr,"showpage\n");
	return; }


/* Execution time is 85 s for itmax=1000, clustered, all bindings; 25 s for
itmax=1000, unclustered, all bindings, 1377 s for itmax=20000, clustered, first
binding. */
void bindrate(void) {
	double x[3],t;
	double tmax=10000000;
	const int itmax=20000;
	const int nhist=200;
	double hist[nhist],delta;
	double lo=-6,hi=7;
	int j,it,bind;
	FILE *fptr;

	delta=(hi-lo)/(nhist-1);
	for(j=0;j<nhist;j++) hist[j]=0;
	for(it=0;it<itmax;) {
		sphererandd(x,1.1*SPHRAD,1.1*SPHRAD);
		bind=diffuse2bind(x,&t);
		if(bind) it++;
		for(t=0;t<tmax&&bind;) {
			move2unbind(x);
			bind=diffuse2bind(x,&t);
			if(bind) {
				bind=0; //*****  Enable this line for 1st binding, remove for all bindings
				j=(int)floor((log10(t)-lo)/delta)+1;
				if(j<0) j=0;
				else if(j>=nhist) j=nhist-1;
				hist[j]+=1.0/itmax; }}}
	fptr=fopen("bindrate","w");
	for(it=0;it<nhist;it++)
		fprintf(fptr,"%lf %lf\n",lo+((double)it-0.5)*delta,hist[it]);
	fclose(fptr);
return; }





void outputrecept(void) {
	FILE *fptr;
	double latmin;

	fptr=NULL;
//	fptr=fopen("latticepts","w");
	if(CLUSTERED) latmin=PI/2.0-CLRAD/SPHRAD;
	else latmin=-PI/2.0;
	spherelatticerows(NRECEPT,SPHRAD,latmin,fptr);
//	fclose(fptr);
	return; }


void teststuff(void) {
	int i;
	double x[3],xr[3];

	x[0]=x[1]=x[2]=0;
	for(i=0;i<1000;i++) {
		printf("enter x vector: ");
		scanf("%lf %lf %lf",&x[0],&x[1],&x[2]);
		xr[0]=xr[1]=xr[2]=0;
		printf("%lf %lf  %lf %lf %lf\n",dist2surf(x),dist2recept(x,xr),xr[0],xr[1],xr[2]); }
	return; }




/* ********************** MAIN SEGMENT *********************** */


int main(void) {
	time_t tstt;

	tstt=time(NULL);
	randomize();
	randtable(GAUSS,4096,1);
	printf("Program starting.\n");

//	outputrecept();
	teststuff();
//	geminatebind();
//	equilibrium();
//	initialsite();
//	allsites();
//	radcorrfn();
//	radcorrall();
//	barfigure();
//	bindrate();
	printf("execution time: %1.0f seconds\n",difftime(time(NULL),tstt));

	return 0; }