%------------------------------------------------------------
% 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];

% Counted everything
% raw_data = [10  579 60.50;
%             15  415 1146.32;
%             20  663 3959.48;
%             30  262 10370;
%             35  61  1394.16;
%             40  336 25278;
%             55  28  10000;
%             60  122 58876;
%             150 126 219493];
        
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),'.'); hold on;
vherrorbar(data(:,1),data(:,2),ones(N,1),data(:,3));
hold on;
%herrorbar(data(:,1),data(:,2),ones(size(data(:,1))),'.');
xlabel('Incident angle (in degrees)');
ylabel('Rate of scattering (counts per second)');
title('Fits to weighted (red) and vanilla (black) Rutherford diff. cross sections');

%------------------------------------------------------------
% FITTING: Rutherford cross section. NOTE that in the fit
%   the angles are converted into radians.
%------------------------------------------------------------
sigx = pi/180*1*zeros(N,1);
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),sigx,...
    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),sigx,...
    @weighted_rdcs,a0,dYda,sgn);
%plot(data(:,1),yfit_rdcs,'k','linewidth',1);
%plot(data(:,1),yfit_w_rdcs,'r','linewidth',2);

% Smooth curve drawn
theta = linspace(10,150,1000);
rdcs = function_rdcs(pi/180*theta,a_rdcs);
wrdcs = weighted_rdcs(pi/180*theta,a_w_rdcs);
hold on;
plot(theta,wrdcs,'r','linewidth',2);
plot(theta,rdcs,'k--');