%------------------------------------------------------------
% Tony Hyun Kim, 9/30/2008
% Data is in the following format: [angle(deg), count, time(s)]
%------------------------------------------------------------
       
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;
hold on;

%------------------------------------------------------------
% FITTING: Rutherford cross section. NOTE that in the fit
%   the angles are converted into radians.
%------------------------------------------------------------
sigx = pi/180*1*ones(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);

ag0 = [100 1];
function_gaussian = @(x,a) a(1)*exp(-a(2)*x.^2);
[a_g,aerr_g,chisq_g,yfit_g,corr] = levmar(pi/180*data(:,1),data(:,2),data(:,3),sigx,...
    function_gaussian,ag0,dYda,sgn);


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

fs = 15;

xlabel('Howitzer angle \phi (deg)','fontsize',fs);
ylabel('Rate of scattering (cps)','fontsize',fs);
title('Angular dependence of the \alpha-particle scattering rate','fontsize',fs);

legend('Measured data','Convolved cross-section (\chi^2/dof = 0.93)',...
        'Original Rutherford cross-section (\chi^2/dof = 12.5)',1);


% Illustrate additional y-errors from the uncertainty in x
x = pi/180*[10 15 20 30 40 55 60 145]';

Nx = length(x);
x_diff = zeros(1,Nx-1);
for i = 1:(Nx-1)
    x_diff(i) = x(i+1)-x(i);
end
dx = min(x_diff)/1000;

ydy_eval = weighted_rdcs(x+dx,a_w_rdcs);
y_eval   = weighted_rdcs(x-dx,a_w_rdcs);
dydx_eval = (ydy_eval-y_eval)/(2*dx);
dydx_abs  = abs(dydx_eval);
sig = data(:,3);
eff_sig = sqrt(sig.^2 + (dydx_abs.*sigx).^2);

vherrorbar(data(:,1),data(:,2),ones(N,1),data(:,3));
%vherrorbar(data(2:Nx,1),data(2:Nx,2),ones(N-1,1),eff_sig(2:Nx));

