ww = new double[ndat];
final double POTN=1.571000,BIG=1.0e30,ACC=1.0e-6;
final double PI=3.141592653589793238;
Gamma gam = new Gamma();
Brent brent = new Brent(ACC);
Chixy chixy = new Chixy();
int j;
double amx,amn,scale,bmn,bmx,d1,d2,r2;
double[] ang = new double[7];
double[] ch = new double[7];
doubleW varx = new doubleW(0);
doubleW vary = new doubleW(0);
doubleW dum1 = new doubleW(0);
Moment.avevar(x,dum1,varx);
Moment.avevar(y,dum1,vary);
scale=sqrt(varx.val/vary.val);
for (j=0;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]));
}
Fitab fit = new Fitab(xx,yy,ww);
b = fit.b;
offs=ang[0]=0.0;
ang[1]=atan(b);
ang[3]=0.0;
ang[4]=ang[1];
ang[5]=POTN;
for (j=3;j<6;j++) ch[j]=chixy.get(ang[j]);
brent.bracket(ang[0],ang[1],chixy);
ang[0] = brent.ax; ang[1] = brent.bx; ang[2] = brent.cx;
ch[0] = brent.fa; ch[1] = brent.fb; ch[2] = brent.fc;
b = brent.minimize(chixy);
chi2=chixy.get(b);
a=aa;
q=gam.gammq(0.5*(ndat-2),chi2*0.5);
r2=0.0;
for (j=0;j<ndat;j++) r2 += ww[j];