#include <math.h>#define NRANSI#include "nrutil.h"#define POTN 1.571000#define BIG 1.0e30#define PI 3.14159265#define ACC 1.0e-3int nn;float *xx,*yy,*sx,*sy,*ww,aa,offs;void fitexy(float x[], float y[], int ndat, float sigx[], float sigy[],	float *a, float *b, float *siga, float *sigb, float *chi2, float *q){	void avevar(float data[], unsigned long n, float *ave, float *var);	float brent(float ax, float bx, float cx,		float (*f)(float), float tol, float *xmin);	float chixy(float bang);	void fit(float x[], float y[], int ndata, float sig[], int mwt,		float *a, float *b, float *siga, float *sigb, float *chi2, float *q);	float gammq(float a, float x);	void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,		float *fc, float (*func)(float));	float zbrent(float (*func)(float), float x1, float x2, float tol);	int j;	float swap,amx,amn,varx,vary,ang[7],ch[7],scale,bmn,bmx,d1,d2,r2,		dum1,dum2,dum3,dum4,dum5;	xx=vector(1,ndat);	yy=vector(1,ndat);	sx=vector(1,ndat);	sy=vector(1,ndat);	ww=vector(1,ndat);	avevar(x,ndat,&dum1,&varx);	avevar(y,ndat,&dum1,&vary);	scale=sqrt(varx/vary);	nn=ndat;	for (j=1;j<=ndat;j++) {		xx[j]=x[j];		yy[j]=y[j]*scale;		sx[j]=sigx[j];		sy[j]=sigy[j]*scale;		ww[j]=sqrt(SQR(sx[j])+SQR(sy[j]));	}	fit(xx,yy,nn,ww,1,&dum1,b,&dum2,&dum3,&dum4,&dum5);	offs=ang[1]=0.0;	ang[2]=atan(*b);	ang[4]=0.0;	ang[5]=ang[2];	ang[6]=POTN;	for (j=4;j<=6;j++) ch[j]=chixy(ang[j]);	mnbrak(&ang[1],&ang[2],&ang[3],&ch[1],&ch[2],&ch[3],chixy);	*chi2=brent(ang[1],ang[2],ang[3],chixy,ACC,b);	*chi2=chixy(*b);	*a=aa;	*q=gammq(0.5*(nn-2),*chi2*0.5);	for (r2=0.0,j=1;j<=nn;j++) r2 += ww[j];	r2=1.0/r2;	bmx=BIG;	bmn=BIG;	offs=(*chi2)+1.0;	for (j=1;j<=6;j++) {		if (ch[j] > offs) {			d1=fabs(ang[j]-(*b));			while (d1 >= PI) d1 -= PI;			d2=PI-d1;			if (ang[j] < *b) {				swap=d1;				d1=d2;				d2=swap;			}			if (d1 < bmx) bmx=d1;			if (d2 < bmn) bmn=d2;		}	}	if (bmx < BIG) {		bmx=zbrent(chixy,*b,*b+bmx,ACC)-(*b);		amx=aa-(*a);		bmn=zbrent(chixy,*b,*b-bmn,ACC)-(*b);		amn=aa-(*a);		*sigb=sqrt(0.5*(bmx*bmx+bmn*bmn))/(scale*SQR(cos(*b)));		*siga=sqrt(0.5*(amx*amx+amn*amn)+r2)/scale;	} else (*sigb)=(*siga)=BIG;	*a /= scale;	*b=tan(*b)/scale;	free_vector(ww,1,ndat);	free_vector(sy,1,ndat);	free_vector(sx,1,ndat);	free_vector(yy,1,ndat);	free_vector(xx,1,ndat);}#undef POTN#undef BIG#undef PI#undef ACC#undef NRANSI