1 Text format
using System;
namespace Legalsoft.Truffer
{
public class Fitexy
{
private double a { get; set; }
private double b { get; set; }
private double siga { get; set; }
private double sigb { get; set; }
private double chi2 { get; set; }
private double q { get; set; }
private int ndat { get; set; }
private double[] xx { get; set; }
private double[] yy { get; set; }
private double[] sx { get; set; }
private double[] sy { get; set; }
private double[] ww { get; set; }
private double aa;
private double offs;
public Fitexy(double[] x, double[] y, double[] sigx, double[] sigy)
{
this.ndat = x.Length;
this.xx = new double[ndat];
this.yy = new double[ndat];
this.sx = new double[ndat];
this.sy = new double[ndat];
this.ww = new double[ndat];
const double POTN = 1.571000;
const double BIG = 1.0e30;
const double ACC = 1.0e-6;
Gamma gam = new Gamma();
//Brent brent = new Brent(ACC);
Chixy chixy = new Chixy(xx, yy, sx, sy, ww, ref aa, ref offs);
double varx = 0.0;
double vary = 0.0;
double[] ang = new double[7];
double[] ch = new double[7];
double dum1 = 0.0;
Moment.avevar(x, ref dum1, ref varx);
Moment.avevar(y, ref dum1, ref vary);
double scale = Math.Sqrt(varx / vary);
for (int 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] = Math.Sqrt(Globals.SQR(sx[j]) + Globals.SQR(sy[j]));
}
Fitab fit = new Fitab(xx, yy, ww);
b = fit.b;
offs = ang[0] = 0.0;
ang[1] = Math.Atan(b);
ang[3] = 0.0;
ang[4] = ang[1];
ang[5] = POTN;
for (int j = 3; j < 6; j + + )
{
ch[j] = chixy.funk(ang[j]);
}
Brent brent = new Brent(ACC);
//brent.func = chixy.functorMethod;
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.funk(b);
a = aa;
q = gam.gammq(0.5 * (ndat – 2), chi2 * 0.5);
double r2 = 0.0;
for (int j = 0; j < ndat; j + + )
{
r2 + = ww[j];
}
r2 = 1.0 / r2;
double bmx = BIG;
double bmn = BIG;
offs = chi2 + 1.0;
for (int j = 0; j < 6; j + + )
{
if (ch[j] > offs)
{
double d1 = Math.Abs(ang[j] – b);
while (d1 >= Math.PI)
{
d1 -= Math.PI;
}
double d2 = Math.PI – d1;
if (ang[j] < b)
{
Globals.SWAP(ref d1, ref d2);
}
if (d1 < bmx)
{
bmx = d1;
}
if (d2 < bmn)
{
bmn = d2;
}
}
}
if (bmx < BIG)
{
bmx = Roots.zbrent(chixy, b, b + bmx, ACC) – b;
double amx = aa – a;
bmn = Roots.zbrent(chixy, b, b – bmn, ACC) – b;
double amn = aa – a;
sigb = Math.Sqrt(0.5 * (bmx * bmx + bmn * bmn)) / (scale * Globals.SQR(Math.Cos(b)));
siga = Math.Sqrt(0.5 * (amx * amx + amn * amn) + r2) / scale;
}
else sigb = siga = BIG;
a /= scale;
b = Math.Tan(b) / scale;
}
}
public class Chixy : UniVarRealValueFun
{
private double[] xx { get; set; }
private double[] yy { get; set; }
private double[] sx { get; set; }
private double[] sy { get; set; }
private double[] ww { get; set; }
private double aa { get; set; }
private double offs { get; set; }
public Chixy(double[] xxx, double[] yyy, double[] ssx, double[] ssy, double[] www, ref double aaa, ref double ooffs)
{
this.xx = xxx;
this.yy = yyy;
this.sx = ssx;
this.sy = ssy;
this.ww = www;
this.aa = aaa;
this.offs = ooffs;
}
public double funk(double bang)
{
return get(bang);
}
public double get(double bang)
{
const double BIG = 1.0e30;
//int j;
int nn = xx.Length;
//double ans;
double avex = 0.0;
double avey = 0.0;
double sumw = 0.0;
//double b;
double b = Math.Tan(bang);
for (int j = 0; j < nn; j + + )
{
ww[j] = Globals.SQR(b * sx[j]) + Globals.SQR(sy[j]);
sumw + = (ww[j] = (ww[j] < 1.0 / BIG ? BIG : 1.0 / ww[j]));
avex + = ww[j] * xx[j];
avey + = ww[j] * yy[j];
}
avex /= sumw;
avey /= sumw;
aa = avey – b * avex;
double ans = -offs;
for (int j = 0; j < nn; j + + )
{
ans + = ww[j] * Globals.SQR(yy[j] – aa – b * xx[j]);
}
return ans;
}
}
}
2 Code format
using System; namespace Legalsoft.Truffer { public class Fitexy { private double a { get; set; } private double b { get; set; } private double siga { get; set; } private double sigb { get; set; } private double chi2 { get; set; } private double q { get; set; } private int ndat { get; set; } private double[] xx { get; set; } private double[] yy { get; set; } private double[] sx { get; set; } private double[] sy { get; set; } private double[] ww { get; set; } private double aa; private double offs; public Fitexy(double[] x, double[] y, double[] sigx, double[] sigy) { this.ndat = x.Length; this.xx = new double[ndat]; this.yy = new double[ndat]; this.sx = new double[ndat]; this.sy = new double[ndat]; this.ww = new double[ndat]; const double POTN = 1.571000; const double BIG = 1.0e30; const double ACC = 1.0e-6; Gamma gam = new Gamma(); //Brent brent = new Brent(ACC); Chixy chixy = new Chixy(xx, yy, sx, sy, ww, ref aa, ref offs); double varx = 0.0; double vary = 0.0; double[] ang = new double[7]; double[] ch = new double[7]; double dum1 = 0.0; Moment.avevar(x, ref dum1, ref varx); Moment.avevar(y, ref dum1, ref vary); double scale = Math.Sqrt(varx / vary); for (int 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] = Math.Sqrt(Globals.SQR(sx[j]) + Globals.SQR(sy[j])); } Fitab fit = new Fitab(xx, yy, ww); b = fit.b; offs = ang[0] = 0.0; ang[1] = Math.Atan(b); ang[3] = 0.0; ang[4] = ang[1]; ang[5] = POTN; for (int j = 3; j < 6; j + + ) { ch[j] = chixy.funk(ang[j]); } Brent brent = new Brent(ACC); //brent.func = chixy.functorMethod; 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.funk(b); a = aa; q = gam.gammq(0.5 * (ndat - 2), chi2 * 0.5); double r2 = 0.0; for (int j = 0; j < ndat; j + + ) { r2 + = ww[j]; } r2 = 1.0 / r2; double bmx = BIG; double bmn = BIG; offs = chi2 + 1.0; for (int j = 0; j < 6; j + + ) { if (ch[j] > offs) { double d1 = Math.Abs(ang[j] - b); while (d1 >= Math.PI) { d1 -= Math.PI; } double d2 = Math.PI - d1; if (ang[j] < b) { Globals.SWAP(ref d1, ref d2); } if (d1 < bmx) { bmx = d1; } if (d2 < bmn) { bmn = d2; } } } if (bmx < BIG) { bmx = Roots.zbrent(chixy, b, b + bmx, ACC) - b; double amx = aa - a; bmn = Roots.zbrent(chixy, b, b - bmn, ACC) - b; double amn = aa - a; sigb = Math.Sqrt(0.5 * (bmx * bmx + bmn * bmn)) / (scale * Globals.SQR(Math.Cos(b))); siga = Math.Sqrt(0.5 * (amx * amx + amn * amn) + r2) / scale; } else sigb = siga = BIG; a /= scale; b = Math.Tan(b) / scale; } } public class Chixy : UniVarRealValueFun { private double[] xx { get; set; } private double[] yy { get; set; } private double[] sx { get; set; } private double[] sy { get; set; } private double[] ww { get; set; } private double aa { get; set; } private double offs { get; set; } public Chixy(double[] xxx, double[] yyy, double[] ssx, double[] ssy, double[] www, ref double aaa, ref double ooffs) { this.xx = xxx; this.yy = yyy; this.sx = ssx; this.sy = ssy; this.ww = www; this.aa = aaa; this.offs = ooffs; } public double funk(double bang) { return get(bang); } public double get(double bang) { const double BIG = 1.0e30; //int j; int nn = xx.Length; //double ans; double avex = 0.0; double avey = 0.0; double sumw = 0.0; //double b; double b = Math.Tan(bang); for (int j = 0; j < nn; j + + ) { ww[j] = Globals.SQR(b * sx[j]) + Globals.SQR(sy[j]); sumw + = (ww[j] = (ww[j] < 1.0 / BIG ? BIG : 1.0 / ww[j])); avex + = ww[j] * xx[j]; avey + = ww[j] * yy[j]; } avex /= sumw; avey /= sumw; aa = avey - b * avex; double ans = -offs; for (int j = 0; j < nn; j + + ) { ans + = ww[j] * Globals.SQR(yy[j] - aa - b * xx[j]); } return ans; } } }