clear all;
close all;
format long;

e = 1.602176634e-19;
U = 35000;
c = 299792458;
h = 6.62607015e-34;
m = 9.1093837015e-31;
E0 = 510.99895000;  %keV, Ruheenergie des Elektrons

data = readmatrix('messung_2020-05-19_Compton.txt');
ch = data(:,1);
N000 = data(:,2);
N030 = data(:,8);
N045 = data(:,10);
N060 = data(:,12);
N075 = data(:,14);
N090 = data(:,16);
N105 = data(:,18);
N120 = data(:,20);
N135 = data(:,22);

%Kalibration
a = 0.07081;
b = -0.5051;
E = a.*ch + b;

figure;
grid.LineWidth = 1;
ax = gca;
ax.FontSize = 10;
ax.LineWidth = 1.5;
plot(E, N000, 'linewidth', 1.5, 'Color', '#EDB120');
hold on;
plot(E, N060, 'linewidth', 1.5, 'Color', [0 169/255 224/255]);
hold on;
plot(E, N135, 'linewidth', 1.5, 'Color', [166/255 225/255 244/255]);
xlim([0 35]);
xlabel('Energie der Photonen, E_\gamma (keV)', 'FontSize', 12);
ylabel('Pulse', 'FontSize', 12);
legend('0°', '60°', '135°', 'FontSize', 10);
legend('boxoff');

%Peakpositionen
%ohne 0deg-Spektrum (weil das schon bei der Kalibration verwendet wurde)
theta = [30 45 60 75 90 105 120 135];
E_Kalpha  = [17.48 17.41 17.27 17.13 16.99 16.70 16.63 16.56];
dE_Kalpha = [ 0.07  0.07  0.07  0.07  0.07  0.14  0.14  0.07];
w = 1./dE_Kalpha.^2;

figure;
grid.LineWidth = 1;
ax = gca;
ax.FontSize = 10;
ax.LineWidth = 1.5;

%theoretischer Verlauf
theta_theo = [0:140];
E_theo = E_Kalpha(1)./( 1 + E_Kalpha(1)/E0.*(1-cosd(theta_theo)) );
l1 = plot(theta_theo, E_theo, 'Color', '#F5D589', 'linewidth', 1.5);
hold on;
p1 = errorbar([0], [17.48], [0.07], '+', 'Color', '#EDB120', 'linewidth', 1.5);
hold on;

%Fit mit Ruheenergie des Elektrons als freier Parameter
ft = fittype( @(E_gamma, E_elektron, x) E_gamma./(1 + E_gamma./E_elektron.*(1-cosd(x))));
[curve, gof] = fit(theta', E_Kalpha', ft, 'Start', [17.48, 511], 'Weights', w)

fit_output = curve.E_gamma./(1+curve.E_gamma./curve.E_elektron.*(1-cosd(theta_theo)));
l2 = plot(theta_theo, fit_output, 'Color', [166/255 225/255 244/255], 'linewidth', 1.5);
hold on;

%Messwerte
p2 = errorbar(theta, E_Kalpha, dE_Kalpha, '+', 'Color', [0 169/255 224/255], 'linewidth', 1.5);
hold on;

xlabel('Streuwinkel, \theta (°)', 'FontSize', 12);
ylabel('Energie, E_\gamma^\prime (keV)', 'FontSize', 12);
%legend([p1 l1], {'E_\gamma bei 0°', 'Theorie'}, 'FontSize', 10);
%legend([p1 l1 p2], {'E_\gamma bei 0°', 'Theorie', 'Messwerte'}, 'FontSize', 10);
legend([p1 l1 p2 l2], {'E_\gamma bei 0°', 'Theorie', 'Messwerte', 'Fit'}, 'FontSize', 10);
legend('boxoff');

%L = fittype( @(C, x0, g, x) C./( (x.^2 - x0^2).^2 + g^2*x0^2))
%L_test = fittype( @(C, x) C./( (x.^2 - 254^2).^2 + 4^2*254^2))
%f = fit( E(249:260), N060(249:260), L, 'Start', [1.5e9, 4, 17.4] )
%plot( f, E(249:260), N060(249:260) )

hold off;