clear all;

%histogramm
%daten einlesen
d=load('histogram.txt');
d=d';
U=596;

%korrektur:
k=[1:numel(d)];
gt(k)=1./(1 - d(k).*(1.6.*10.^(-004)));
d=d.*gt;

cps = [1050:10:1300]';
cps_n = zeros(25,1);

for i=[1:25];
    for k=[1:numel(d)]
        if (d(k) > cps(i)) & (d(k) <= cps(i+1))
            cps_n(i) = cps_n(i)+1;
        end
    end
end
i=[1:25];

cps_n_f = sqrt(cps_n);
cps_n_f(cps_n_f==0) = 1;
cps_d = cps+5;              %um halbe balkenbreite korrigiert

%gausfit
% n=a(1) * exp(-0.5* (cps-a(2)).^2/a(3)^2);
    a=[35,1000,30];
	eps = 0.05;     %max relativer fehler
	
	ndata = numel(i);	
	ma = 3;
	
	[a,sd_a,var,covar,c] = mrqmin(cps_d(i),cps_n(i),cps_n_f(i),ndata,ma,a,eps,'funcs');
    disp('n=a(1) * exp(-0.5* (cps-a(2)).^2/a(3)^2)');
    disp(a);
    disp(sd_a);
    

    cps_x = [min(cps_d(i)):1:max(cps_d(i))];
    n=a(1) * exp(-0.5* (cps_x-a(2)).^2/a(3)^2);

    
f1=figure(1);clf;
hold on;
p1=bar(cps_d(i),cps_n(i),'r');
errorbar(cps_d(i),cps_n(i),cps_n_f(i),'b.');
p3=plot(cps_x,n);

ax1=gca;
%set(ax1,'YLim',[0,max(Pld)*1.1]);
%set(ax1,'XLim',[0,max(Id)*1.1]); e^{- (z-z_0)^2/ \sigma^2}}

legend([p1,p3],'Gemessen N','Fit N=a \cdot e^{- (z-z_0)^2/ \sigma^2}',2);
xlabel('z');
ylabel('N');

hold off;

print -deps2 -tiff -r300 histogram


for j=[1:8]
    for i=[1:37]
        if j*i <= numel(d)
            dd(i,j)=d(j*i);
        else
            dd(i,j)=NaN;
        end
    end
end

        

t=['$z_{unkorigiert}$ & & & & & & &\\'];
d=latex(dd,'%3.1f','nomath');
f = fopen('histogram.tab','w');
fwrite(f,d);
fclose(f);
f = fopen('histogram_h.tab','w');
fwrite(f,t);
fclose(f);

