%------------------------------------------------------------
% Tony Hyun Kim, 9/30/2008
% Data is in the following format: [angle(deg), count, time(s)]
%------------------------------------------------------------
raw_data = [ 10     377     60.50;
             20     395     3959.48;
             30     146     10370;
             40      18     5278.2];

N = length(raw_data(:,1));
data = zeros(N,3); % [angle(deg), rate, error]

for i = 1:N
    data(i,1) = raw_data(i,1);
    data(i,2) = raw_data(i,2)/raw_data(i,3);
    data(i,3) = sqrt(raw_data(i,2))/raw_data(i,3);
end

figure;
semilogy(data(:,1),data(:,2),'o','markersize',4); hold on;
errorbar(data(:,1),data(:,2),data(:,3),'o','markersize',4);
xlabel('Incident angle');
ylabel('Rate of scattering');

%------------------------------------------------------------
% FITTING: Rutherford cross section. NOTE that in the fit
%   the angles are converted into radians.
%------------------------------------------------------------
function_rdcs = @(x,a) a(1)*(sin(x/2)).^(-4);
a0 = [sin(pi/180*10)^4*data(1,2)];
dYda = {@(x,a) 1}; sgn = 0;
[a_rdcs,aerr,chisq,yfit_rdcs,corr] = levmar(pi/180*data(:,1),data(:,2),data(:,3),...
    function_rdcs,a0,dYda,sgn);
[a_w_rdcs,aerr_w,chisq_w,yfit_w_rdcs,corr] = levmar(pi/180*data(:,1),data(:,2),data(:,3),...
    @weighted_rdcs,a0,dYda,sgn);
plot(data(:,1),yfit_rdcs,'k','linewidth',1);
plot(data(:,1),yfit_w_rdcs,'r','linewidth',2);
title('Fit to weighted (red) and vanilla (green) Rutherford diff. cross section');