//---------------------------------------------------------------------------
#pragma hdrstop

#include "Recipes.h"

float *circx,*circy,*circw;
int circn;

//---------------------------------------------------------------------------
#pragma package(smart_init)
#include <math.h>
#define NRANSI#define EPS 1.0e-7void medfit(float x[], float y[], int ndata, float *a, float *b, float *abdev){	float rofunc(float b);	int j;	float bb,b1,b2,del,f,f1,f2,sigb,temp;	float sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,chisq=0.0;	ndatat=ndata;	xt=x;	yt=y;	for (j=1;j<=ndata;j++) {		sx += x[j];		sy += y[j];		sxy += x[j]*y[j];		sxx += x[j]*x[j];	}	del=ndata*sxx-sx*sx;	aa=(sxx*sy-sx*sxy)/del;	bb=(ndata*sxy-sx*sy)/del;	for (j=1;j<=ndata;j++)		chisq += (temp=y[j]-(aa+bb*x[j]),temp*temp);	sigb=sqrt(chisq/del);	b1=bb;	f1=rofunc(b1);	b2=bb+SIGN(3.0*sigb,f1);	f2=rofunc(b2);	if (b2 == b1) {		*a=aa;		*b=bb;		*abdev=abdevt/ndata;		return;	}	while (f1*f2 > 0.0) {		bb=b2+1.6*(b2-b1);		b1=b2;		f1=f2;		b2=bb;		f2=rofunc(b2);	}	sigb=0.01*sigb;	while (fabs(b2-b1) > sigb) {		bb=b1+0.5*(b2-b1);		if (bb == b1 || bb == b2) break;		f=rofunc(bb);		if (f*f1 >= 0.0) {			f1=f;			b1=bb;		} else {			f2=f;			b2=bb;		}	}	*a=aa;	*b=bb;	*abdev=abdevt/ndata;}
float rofunc(float b){	//float select(unsigned long k, unsigned long n, float arr[]);	int j;	float *arr,d,sum=0.0;	arr = new float[ndatat+1];	for (j=1;j<=ndatat;j++)        {            arr[j]=yt[j]-b*xt[j];        }	if (ndatat & 1) {		aa=select((ndatat+1)>>1,ndatat,arr);	}	else {		j=ndatat >> 1;		aa=0.5*(select(j,ndatat,arr)+select(j+1,ndatat,arr));	}	abdevt=0.0;	for (j=1;j<=ndatat;j++) {		d=yt[j]-(b*xt[j]+aa);		abdevt += fabs(d);		if (yt[j] != 0.0) d /= fabs(yt[j]);		if (fabs(d) > EPS) sum += (d >= 0.0 ? xt[j] : -xt[j]);	}        delete arr;	return sum;}float select(unsigned long k, unsigned long n, float arr[]){	unsigned long i,ir,j,l,mid;	float a,temp;	l=1;	ir=n;	for (;;) {		if (ir <= l+1) {			if (ir == l+1 && arr[ir] < arr[l]) {				SWAP(arr[l],arr[ir])			}			return arr[k];		} else {			mid=(l+ir) >> 1;			SWAP(arr[mid],arr[l+1])			if (arr[l] > arr[ir]) {				SWAP(arr[l],arr[ir])			}			if (arr[l+1] > arr[ir]) {				SWAP(arr[l+1],arr[ir])			}			if (arr[l] > arr[l+1]) {				SWAP(arr[l],arr[l+1])			}			i=l+1;			j=ir;			a=arr[l+1];			for (;;) {				do i++; while (arr[i] < a);				do j--; while (arr[j] > a);				if (j < i) break;				SWAP(arr[i],arr[j])			}			arr[l+1]=arr[j];			arr[j]=a;			if (j >= k) ir=j-1;			if (j <= k) l=i;		}	}}#undef SWAP#undef EPS#undef NRANSI


float fcircfit(float p[])
{
        int i;
        float val = 0;
        float val1,val2;

        for (i=1;i<=circn;i++)
        {
                val1 = sqrt(pow(circx[i]-p[1],2)+pow(circy[i]-p[2],2));
                val2 = val1 - p[3];
                val += circw[i] * pow(val2,2);
        }

        return val;
}

void fcircderiv(float p[], float xi[])
{
        int i;
        float val1,val2;

        xi[1] = 0;
        xi[2] = 0;
        xi[3] = 0;

        for (i=1;i<=circn;i++)
        {
                val1 = sqrt(pow(circx[i]-p[1],2)+pow(circy[i]-p[2],2));
                val2 = val1 - p[3];

                if (val1 > 0)
                {

                xi[1] += circw[i] * 2*p[1]*val2*(circx[i] - p[1])/val1;
                xi[2] += circw[i] * 2*p[2]*val2*(circy[i] - p[2])/val1;
                xi[3] += circw[i] * 2 * val2;
                }
        }
}

#define ITMAX 200
#define EPS 1.0e-10#define FREEALL free_vector(xi,1,n);free_vector(h,1,n);free_vector(g,1,n);void frprmn(float p[], int n, float ftol, int *iter, float *fret,	float (*func)(float []), void (*dfunc)(float [], float [])){	void linmin(float p[], float xi[], int n, float *fret,		float (*func)(float []));	int j,its;	float gg,gam,fp,dgg;	float *g,*h,*xi;	g=vector(1,n);	h=vector(1,n);	xi=vector(1,n);	fp=(*func)(p);	(*dfunc)(p,xi);	for (j=1;j<=n;j++) {		g[j] = -xi[j];		xi[j]=h[j]=g[j];	}	for (its=1;its<=ITMAX;its++) {		*iter=its;		linmin(p,xi,n,fret,func);		if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {			FREEALL			return;		}		fp=(*func)(p);		(*dfunc)(p,xi);		dgg=gg=0.0;		for (j=1;j<=n;j++) {			gg += g[j]*g[j];			dgg += (xi[j]+g[j])*xi[j];		}		if (gg == 0.0) {			FREEALL			return;		}		gam=dgg/gg;		for (j=1;j<=n;j++) {			g[j] = -xi[j];			xi[j]=h[j]=g[j]+gam*h[j];		}	}	nrerror("Too many iterations in frprmn");}#undef ITMAX#undef EPS#undef FREEALL#define TOL 2.0e-4void linmin(float p[], float xi[], int n, float *fret, float (*func)(float [])){	float brent(float ax, float bx, float cx,		float (*f)(float), float tol, float *xmin);	float f1dim(float x);	void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,		float *fc, float (*func)(float));	int j;	float xx,xmin,fx,fb,fa,bx,ax;	ncom=n;	pcom=vector(1,n);	xicom=vector(1,n);	nrfunc=func;	for (j=1;j<=n;j++) {		pcom[j]=p[j];		xicom[j]=xi[j];	}	ax=0.0;	xx=1.0;	mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);	*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);	for (j=1;j<=n;j++) {		xi[j] *= xmin;		p[j] += xi[j];	}	free_vector(xicom,1,n);	free_vector(pcom,1,n);}

#undef TOL
#define ITMAX 100
#define CGOLD 0.3819660#define ZEPS 1.0e-10#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);float brent(float ax, float bx, float cx, float (*f)(float), float tol,	float *xmin){	int iter;	float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;	float e=0.0;	a=(ax < cx ? ax : cx);	b=(ax > cx ? ax : cx);	x=w=v=bx;	fw=fv=fx=(*f)(x);	for (iter=1;iter<=ITMAX;iter++) {		xm=0.5*(a+b);		tol2=2.0*(tol1=tol*fabs(x)+ZEPS);		if (fabs(x-xm) <= (tol2-0.5*(b-a))) {			*xmin=x;			return fx;		}		if (fabs(e) > tol1) {			r=(x-w)*(fx-fv);			q=(x-v)*(fx-fw);			p=(x-v)*q-(x-w)*r;			q=2.0*(q-r);			if (q > 0.0) p = -p;			q=fabs(q);			etemp=e;			e=d;			if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))				d=CGOLD*(e=(x >= xm ? a-x : b-x));			else {				d=p/q;				u=x+d;				if (u-a < tol2 || b-u < tol2)					d=SIGN(tol1,xm-x);			}		} else {			d=CGOLD*(e=(x >= xm ? a-x : b-x));		}		u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));		fu=(*f)(u);		if (fu <= fx) {			if (u >= x) a=x; else b=x;			SHFT(v,w,x,u)			SHFT(fv,fw,fx,fu)		} else {			if (u < x) a=u; else b=u;			if (fu <= fw || w == x) {				v=w;				w=u;				fv=fw;				fw=fu;			} else if (fu <= fv || v == x || v == w) {				v=u;				fv=fu;			}		}	}	nrerror("Too many iterations in brent");}

#undef ITMAX
#undef CGOLD
#undef ZEPS
#undef SHFT

#define GOLD 1.618034
#define GLIMIT 100.0#define TINY 1.0e-20#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc,	float (*func)(float)){	float ulim,u,r,q,fu,dum;	*fa=(*func)(*ax);	*fb=(*func)(*bx);	if (*fb > *fa) {		SHFT(dum,*ax,*bx,dum)		SHFT(dum,*fb,*fa,dum)	}	*cx=(*bx)+GOLD*(*bx-*ax);	*fc=(*func)(*cx);	while (*fb > *fc) {		r=(*bx-*ax)*(*fb-*fc);		q=(*bx-*cx)*(*fb-*fa);		u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/			(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));		ulim=(*bx)+GLIMIT*(*cx-*bx);		if ((*bx-u)*(u-*cx) > 0.0) {			fu=(*func)(u);			if (fu < *fc) {				*ax=(*bx);				*bx=u;				*fa=(*fb);				*fb=fu;				return;			} else if (fu > *fb) {				*cx=u;				*fc=fu;				return;			}			u=(*cx)+GOLD*(*cx-*bx);			fu=(*func)(u);		} else if ((*cx-u)*(u-ulim) > 0.0) {			fu=(*func)(u);			if (fu < *fc) {				SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))				SHFT(*fb,*fc,fu,(*func)(u))			}		} else if ((u-ulim)*(ulim-*cx) >= 0.0) {			u=ulim;			fu=(*func)(u);		} else {			u=(*cx)+GOLD*(*cx-*bx);			fu=(*func)(u);		}		SHFT(*ax,*bx,*cx,u)		SHFT(*fa,*fb,*fc,fu)	}}#undef GOLD#undef GLIMIT#undef TINY#undef SHFT

float f1dim(float x)
{	int j;	float f,*xt;	xt=vector(1,ncom);	for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];	f=(*nrfunc)(xt);	free_vector(xt,1,ncom);	return f;}

void kstwo(float data1[], unsigned long n1, float data2[], unsigned long n2,
	float *d, float *prob){	float probks(float alam);	void sort(unsigned long n, float arr[]);	unsigned long j1=1,j2=1;	float d1,d2,dt,en1,en2,en,fn1=0.0,fn2=0.0;	sort(n1,data1);	sort(n2,data2);	en1=n1;	en2=n2;	*d=0.0;	while (j1 <= n1 && j2 <= n2) {		if ((d1=data1[j1]) <= (d2=data2[j2])) fn1=j1++/en1;		if (d2 <= d1) fn2=j2++/en2;		if ((dt=fabs(fn2-fn1)) > *d) *d=dt;	}	en=sqrt(en1*en2/(en1+en2));	*prob=probks((en+0.12+0.11/en)*(*d));}

#define EPS1 0.001
#define EPS2 1.0e-8float probks(float alam){	int j;	float a2,fac=2.0,sum=0.0,term,termbf=0.0;	a2 = -2.0*alam*alam;	for (j=1;j<=100;j++) {		term=fac*exp(a2*j*j);		sum += term;		if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum) return sum;		fac = -fac;		termbf=fabs(term);	}	return 1.0;}#undef EPS1#undef EPS2

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7#define NSTACK 50void sort(unsigned long n, float arr[]){	unsigned long i,ir=n,j,k,l=1;	int jstack=0,*istack;	float a,temp;	istack=ivector(1,NSTACK);	for (;;) {		if (ir-l < M) {			for (j=l+1;j<=ir;j++) {				a=arr[j];				for (i=j-1;i>=l;i--) {					if (arr[i] <= a) break;					arr[i+1]=arr[i];				}				arr[i+1]=a;			}			if (jstack == 0) break;			ir=istack[jstack--];			l=istack[jstack--];		} else {			k=(l+ir) >> 1;			SWAP(arr[k],arr[l+1])			if (arr[l] > arr[ir]) {				SWAP(arr[l],arr[ir])			}			if (arr[l+1] > arr[ir]) {				SWAP(arr[l+1],arr[ir])			}			if (arr[l] > arr[l+1]) {				SWAP(arr[l],arr[l+1])			}			i=l+1;			j=ir;			a=arr[l+1];			for (;;) {				do i++; while (arr[i] < a);				do j--; while (arr[j] > a);				if (j < i) break;				SWAP(arr[i],arr[j]);			}			arr[l+1]=arr[j];			arr[j]=a;			jstack += 2;			if (jstack > NSTACK) nrerror("NSTACK too small in sort.");			if (ir-i+1 >= j-l) {				istack[jstack]=ir;				istack[jstack-1]=i;				ir=j-1;			} else {				istack[jstack]=j-1;				istack[jstack-1]=l;				l=i;			}		}	}	free_ivector(istack,1,NSTACK);}#undef M#undef NSTACK#undef SWAP#undef NRANSI

#define NRANSI
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;#define M 7#define NSTACK 50

void sort2fi(unsigned long n, float arr[], int brr[])
{	unsigned long i,ir=n,j,k,l=1;	int *istack,jstack=0;	float a,temp;        int b;	istack=ivector(1,NSTACK);	for (;;) {		if (ir-l < M) {			for (j=l+1;j<=ir;j++) {				a=arr[j];				b=brr[j];				for (i=j-1;i>=l;i--) {					if (arr[i] <= a) break;