raw_data = load('hydrogen.txt');
N = length(raw_data(:,1));

data = zeros(N,3);
for i = 1:N
    data(i,1) = 10-i;
    data(i,2) = raw_data(i,2); % Mechanical readings
    data(i,3) = syserr; % Dominated by systematic error
end

% Consider only the data taken with 3600gpmm to simplify analysis
N = N-1;
data = data(1:N,:);


% Need to convert the mechanical readings to wavelengths using Hg data
%   inverting quad. formula ax^2 + bx + c
a = afit(1);
b = afit(2);
c = afit(3);

for i = 1:N
    c_prime = c - data(i,2);
    data(i,2) = (-b + sqrt(b^2-4*a*c_prime))/(2*a);
    
    % Error propagation
    dl_da = (-4*a*c_prime/sqrt(b^2-4*a*c_prime) - 2*(-b+sqrt(b^2-4*a*c_prime)))/(4*a^2);
    dl_db = (-1 + (b^2-4*a*c_prime)^(-1/2)*b)/(2*a);
    dl_dc = -2*a/sqrt(b^2-4*a*c_prime);
    g = [dl_da dl_db dl_dc]';
    
    dl_dmech = -dl_dc;  
    mecherror = dl_dmech^2*syserr^2;
       
    data(i,3) = g'*corr_ruler*g + mecherror;
    data(i,3) = sqrt(data(i,3));
end

figure;
errorbar(data(:,1),data(:,2),data(:,3),'ro','markersize',4); hold on;

% Attempted fit
function_ryd = @(x,a) 1./(a(1)*(1/4 - 1./x.^2));
a0 = [1.096e-3]; % actual Rydberg
dYda = {@(x,a) 1}; sgn = 0;
[A,Aerr,chisq,yfit,corr] = levmar(data(:,1),data(:,2),data(:,3),...
    function_ryd,a0,dYda,sgn);

M = 1000;
x = linspace(data(1,1),data(N,1),M);
y = function_ryd(x,A);
plot(x,y,'b','linewidth',1);

% Results
chisqdof = chisq/5;
fprintf(1,'Normalized chi sq is %f\n', chisqdof);
