%% pwfitexy.m is my MATLAB version of the fitexy program from Numerical Recipes
%% (see p662 Section 15.3) that determines a linear fit to data when there are 
%% errors in both the x and y data.

function [c,ce,chi2,q]=pwfitexy(xv,xve,yv,yve)

% xv = x-vector data.
% xve = x-vector error.
% yv = y-vector data.
% yve = y-vector error.
%
% c = coefficents of fit [y=c(1)*x+c(2)].
% ce = coefficents of errors of fit.
% chi2 = chi square of fit.
% q = probability of chi square fit.

%Make sure all input are columns
xv=xv(:);
yv=yv(:);
xve=xve(:);
yve=yve(:);

%Define required constants
potn=1.571;
big=1e30;
acc=1e-3;

%Declare required GLOBAL VARIABLES
global xx yy sx sy ww aa offs nop

%Number of points
nop=length(xv);

%Determine x and y variances and scale data as required
varx=std(xv)^2;               %x variance
vary=std(yv)^2;               %y variance
scale=sqrt(varx/vary);        %scale factor
xx=xv;                        %scaled x-vector
yy=yv*scale;                  %scaled y-vector
sx=xve;                       %scaled x-vector error
sy=yve*scale;                 %scaled y-vector error
ww=sqrt(sx.^2+sy.^2);         %weight for first fit

%Determine initial estimate for the slope
ci=pwpolfit(xx,yy,ww,1);
b=ci(1);
offs=0;

%Construct several angles for reference points
ang=[0 atan(b) 0 0 atan(b) potn];
ch=zeros(1,6);
for j=4:6
 ch(j)=pwchixy(ang(j));
end

%Locate the chi square minimum
[ang(1),ang(2),ang(3),ch(1),ch(2),ch(3)]=pwmnbrak('pwchixy',ang(1),ang(2));
b=fminbnd('pwchixy',ang(1),ang(3),optimset('tolx',acc,'display','off'));
chi2=pwchixy(b);

%Compute chi square probability
a=aa;
q=gammainc(0.5*chi2,0.5*(nop-2));

%Save the inverse sum of weights at the minimum
r2=1./sum(ww);

%Find standard errors for b as points where chi square equals one
bmx=big;
bmn=big;
offs=chi2+1;
for j=1:6
 if ch(j)>offs
  d1=mod(abs(ang(j)-b),pi);
  d2=pi-d1;
  if ang(j)<b
   swap=d1;
   d1=d2;
   d2=swap;
  end
  if d1<bmx
   bmx=d1;
  end
  if d2<bmn
   bmn=d2;
  end
 end
end

%Determine roots
if bmx<big
 bmx=fzero('pwchixy',[b b+bmx],optimset('tolx',acc,'display','off'))-b;
 amx=aa-a;
 bmn=fzero('pwchixy',[b b-bmn],optimset('tolx',acc,'display','off'))-b;
 amn=aa-a;
 sigb=sqrt(0.5*(bmx^2+bmn^2))/(scale*cos(b)^2);
 siga=sqrt(0.5*(amx^2+amn^2)+r2)/scale;
else
 sigb=big;
 siga=big;
end

%Unscale the answers
a=a/scale;
b=tan(b)/scale;

%Required form of the output
c=[b a];
ce=[sigb siga]; 


