function solution

%
% This is the raw table given in the problem Set
%
data = [
 0.00      1  11.25      3  22.50     12  33.75     26  45.00     46 
 0.75     90  12.00    124  23.25    168  34.50    210  45.75    252 
 1.50    372  12.75    435  24.00    475  35.25    544  46.50    595 
 2.25    624  13.50    746  24.75    908  36.00    929  47.25    962 
 3.00   1115  14.25   1212  25.50   1221  36.75   1392  48.00   1414 
 3.75   1548  15.00   1518  26.25   1642  37.50   1683  48.75   1801 
 4.50   1848  15.75   1808  27.00   1878  38.25   1906  49.50   1975 
 5.25   1983  16.50   2054  27.75   1981  39.00   1931  50.25   1964 
 6.00   1943  17.25   1917  28.50   1865  39.75   1856  51.00   1855 
 6.75   1763  18.00   1783  29.25   1663  40.50   1580  51.75   1514 
 7.50   1434  18.75   1410  30.00   1264  41.25   1272  52.50   1134 
 8.25   1096  19.50   1015  30.75    890  42.00    842  53.25    705 
 9.00    653  20.25    560  31.50    535  42.75    395  54.00    360 
 9.75    267  21.00    232  32.25    184  43.50    125  54.75    118 
10.50     58  21.75     32  33.00     14  44.25      8  55.50     1
];


%
% Convert it to a two-column matrix
%
data = [data(:,1:2) ; data(:,3:4) ; data(:,5:6) ; data(:,7:8) ; data(:,9:10)];


%
% Get x and N vectors from the raw data
%
x = data(:,1);
N = data(:,2);


%
% calculate y and its error assuming that N is the number of counts,
% thus its error is given by sqrt(N)
%
y = 5*N;
sy = 5*sqrt(N);


%
% Make a plot of y vs x with error bars
%
figure(1)
errorbar(x,y,sy,'o')
xlabel('x')
ylabel('y')


%
% Create an array containing the index number of maxima positions
%
n = [8 23 38 53 67];


%
% Print header for the table to be written
%
fprintf(1,'\n');
fprintf(1, 'i        x_i  +-  sx_i      (chi2_dof)\n');
fprintf(1, '--------------------------------------\n');


%
% Loop over the peaks
%
for peak=1:length(n)
  
  %
  % Take two point on the left, two points on the right
  % and create xx, yy, and sigma_yy (total of 5 points)
  % We will fit these 5 points into a 2nd-deg polynomial
  %
  xx  =  x(n(peak)-2:n(peak)+2);
  yy  =  y(n(peak)-2:n(peak)+2);
  syy = sy(n(peak)-2:n(peak)+2);
  
  
  %
  % Create the matrix-M and vector-B as derived in the lecture
  %
  B = zeros(3,1);
  for i=1:3
    for j=1:3
      M(i,j) = sum(xx.^(i+j-2) ./ syy.^2);
    end
    
    B(i) = sum(yy.*xx.^(i-1) ./ syy.^2);
  end
  
  
  %
  % Evaluate inverse of M
  %
  C = inv(M);
  

  %
  % Evaluate the polynomial coefficients
  %
  a = C*B;
  

  %
  % Evaluate the chi2 per d.o.f. for the fit
  %
  chi2dof = sum(((yy - (a(1) + a(2)*xx + a(3)*xx.^2)) ./ syy) .^2) / 2;

  
  %
  % Create a smooth curve using the polynomail coefs
  % staring at first XX, ending at last xx.
  % Super impose it on the plot (hold/on-off -> superimposes)
  fx = xx(1):0.01:xx(5);
  fy = a(1) + a(2)*fx + a(3)*fx.^2;
  
  hold on
  plot(fx,fy,'r') % 'r' is for red curve
  hold off


  %
  % Calculate x_max. (Take derivative, set it to zero, and solve for x_max)
  %
  xmax(peak) = - a(2) / 2 / a(3);

  
  %
  % Calculate the error in x_max.
  % Note that the third term is for coveriance term
  % as a1 and a2 are not completely independent
  %
  sxmax(peak) = sqrt(C(2,2)*(1/(2*a(3)))^2 + C(3,3)*(a(2)/(2*a(3)^2))^2 + C(2,3)*2*abs(a(2)/(4*a(3)^3)));

  
  %
  % Print out the result to be written to the answer sheet
  %
  fprintf(1, '%d    %7.4f  +-  %6.4f     (%6.4f)\n', peak, xmax(peak), sxmax(peak), chi2dof)
    
end


%
% Make an error-bar plot showing the x_max vs peak-number
%
figure(2)
i=1:5;
errorbar(i,xmax, sxmax, 'o')
xlabel('peak')
ylabel('x_{max}')


%
% Find the intercept, slope, their errors, and chi2 per d.o.f
% for a line fit of x_max vs peak-number data
%
S   = sum(1 ./ sxmax.^2);
Sx  = sum(i ./ sxmax.^2);
Sy  = sum(xmax ./ sxmax.^2);
Sxx = sum(i.^2 ./ sxmax.^2);
Sxy = sum(i.*xmax ./ sxmax.^2);

delta = S*Sxx - Sx^2;

n = (Sxx*Sy - Sx*Sxy) / delta;
m = (S*Sxy - Sx*Sy) / delta;

sn = sqrt(Sxx / delta);
sm = sqrt(S / delta);

chi2dof = sum(((xmax - (n + m*i)) ./ sxmax) .^2) / 3;



%
% Make a smooth plot using the fit parameters.
% Superimpose on the errorbar plot
%
fx=1:5;
fy = n + m*fx;

hold on
plot(fx,fy,'r');
hold off


%
% Print out the result
%
fprintf(1,'\n');
fprintf(1,'n        = %.3f +- %.3f\n', n, sn);
fprintf(1,'m        = %.3f +- %.3f\n', m, sm);
fprintf(1,'chi2_dof = %.3f\n', chi2dof);


%
% Note that the period is twice the slope
% due to the square operation in the model (see the Problem Set)
% Find the period and its error accordingly, and print out the result
%
T = 2*m;
sigma_T = 2*sm;

fprintf(1,'\n');
fprintf(1,'T = 2*m = %.3f +- %.3f\n', T, sigma_T);


%
% Find the angular frequency and its error using the period T
%
omega = 2*pi/T;
sigma_omega = 2*pi/T^2 * sigma_T;

fprintf(1,'\n');
fprintf(1,'omega = %.4f +- %.4f\n', omega, sigma_omega);


%
% The data was produced using omega = 0.2827.
% Note that the answer is within one-sigma error.
%
fprintf(1,'\nNote: Real omega = 0.2827\n');


%
% Create figure files in PDF and EPS to be used in LaTeX/Word
%
figure(1)
print -dpdf peaks.pdf
print -depsc2 peaks.eps

figure(2)
print -dpdf xmax.pdf
print -depsc2 xmax.eps

