clear all;

%spektrum
%daten einlesen
%t=10sec
r=0.034;
dr=0.006;
%  m0=9.109 * 10^(-31);  %kg
%  c= 300000000;
m0=0.511;
c=1;


d=load('spektrum1.csv');
B1=d(:,1);
n1=d(:,2:6);

%korrektur:
gt=1./(1 - n1.*(1.6.*10.^(-004)));
n1=n1.*gt*10;

nn1 = 1/5* (n1(:,1)+n1(:,2)+n1(:,3)+n1(:,4)+n1(:,5));
dn1=sqrt(1/5*((nn1-n1(:,1)).^2 + (nn1-n1(:,2)).^2 + (nn1-n1(:,3)).^2 + (nn1-n1(:,4)).^2 + (nn1-n1(:,5)).^2));

d=load('spektrum2.csv');
B2=d(:,1);
n2=d(:,2:11);

%korrektur:
gt=1./(1 - n2.*(1.6.*10.^(-004)));
n2=n2.*gt*10;

nn2 = 1/10* (n2(:,1)+n2(:,2)+n2(:,3)+n2(:,4)+n2(:,5)+n2(:,6)+n2(:,7)+n2(:,8)+n2(:,9)+n2(:,10));
dn2=sqrt(1/10*((nn2-n2(:,1)).^2 + (nn2-n2(:,2)).^2 + (nn2-n2(:,3)).^2 + (nn2-n2(:,4)).^2 + (nn2-n2(:,5)).^2 + (nn2-n2(:,6)).^2 + (nn2-n2(:,7)).^2 + (nn2-n2(:,8)).^2 + (nn2-n2(:,9)).^2 + (nn2-n2(:,10)).^2));

B=[B1;B2];
dB=B.*0.02;
n=[nn1;nn2];
dn=[dn1;dn2];

%hintergrund abziehen
	t=3600+12*60+31;
	zk=40593/t;
	gt=1./(1 - zk.*(1.6.*10.^(-004)));
	nh=zk*gt*10;
	dnh=sqrt(zk*gt*t)/t*10;

	n=n-nh;
	dn = sqrt(dn.^2+dnh.^2);

%impuls
p= 0.3 * r * B;
dp=0.3*r*dB + 0.3*B*dr;
%E = sqrt(p.^2 * c^2 + m0^2*c^4) - m0 * c^2;
E = sqrt(p.^2 * c^2 + m0^2*c^4);
dE = 1/2 .* (p.^2 * c^2 + m0^2*c^4).^(-1/2) .* 2 .* p .* dp;

kurie=sqrt(n)./p.^(3/2);
dkurie = 1/2* n.^(-1/2) ./p.^(3/2) .*dn + 3/2 * sqrt(n)./p.^(5/2) .*dp;


%spektrumfit
% % f=Ka * (Ea0 - E) +  Kb * (Eb0 - E)
%   f=a(1) * (a(2) - E) +  a(3) * (a(4) - E)
     a=[20,1,2,3];
 	eps = 0.000001;     %max relativer fehler
	
	ndata = numel(p);	
	ma = 4;	
	[a,sd_a,var,covar,c] = mrqmin(E,kurie,dkurie,ndata,ma,a,eps,'funcs_spektrum');
     disp('f=a(1) * (a(2) - E) +  a(3) * (a(4) - E)');
     disp(a);
     disp(sd_a);
     disp(var);
     disp(covar);
    
    e = linspace(0,a(4));
    f = a(1) .* (a(2) - e) .* (a(2) > e) +  a(3) .* (a(4) - e) .* (a(4) > e);
    chi2 = sum(1./dkurie.^2 .* (kurie - (a(1) .* (a(2) - E) .* (a(2) > E) +  a(3) .* (a(4) - E) .* (a(4) > E))).^2 );
    disp('chi2');
    disp(chi2);
    

f1=figure(1);clf;
hold on;
p1=errorbar(E,kurie,dkurie,'.');
%herrorbar(E,kurie,dE,'.');
p3=plot(e,f);

grid on;    
ax1=gca;
set(ax1,'YLim',[0,max(kurie)*1.1]);
set(ax1,'XLim',[0,max(e)*1.1]);

le1=legend([p1(1),p3],'n^{1/2} / p^{3/2}','Fit f=K_a \cdot (E_{a0} - E) +  K_b \cdot (E_{b0} - E)',2);
xl1=xlabel('Gesamtenergie E / [MeV]');
yl1=ylabel('n^{1/2} / p^{3/2}');

    
hold off;
print -deps2 -tiff -r300 kurie


f2=figure(2);
hold on;
p1=errorbar(E,n,dn,'.');

grid on;    
ax1=gca;
set(ax1,'YLim',[0,max(n)*1.1]);
set(ax1,'XLim',[0,max(E)*1.1]);

le1=legend(p1(1),'z',2);
xl1=xlabel('Gesamtenergie E / [MeV]');
yl1=ylabel('z');

    
hold off;
print -deps2 -tiff -r300 spektrum







dd=zeros(numel(B1)+numel(B2),10)*NaN;
dd(1:numel(B1),1:5)=n1;
dd(numel(B1)+1:numel(B1)+numel(B2),:)=n2;

dc=[B,dB,n,dn,p,dp,E,dE,kurie,dkurie];

  t=['$B/[mT]$ & $n1$ & $n2$ & $n3$ & $n4$ & $n5$ & $n6$ & $n7$ & $n8$ & $n9$ & $n10$ \\'];
  d=latex([B,dd],'%3.1f','nomath');
  f = fopen('spektrum1.tab','w');
  fwrite(f,d);
  fclose(f);
  f = fopen('spektrum1_h.tab','w');
  fwrite(f,t);
  fclose(f);

  t=['$B/[mT]$ & $\Delta B$ & $n$ & $\Delta n$ & $p/[MeV]$ & $\Delta p$ & $E/[MeV]$ & $\Delta E$ & $k/[MeV]^{\frac{3}{2}}$ & $\Delta k$\\'];
  d=latex(dc,'%3.1f','nomath');
  f = fopen('spektrum2.tab','w');
  fwrite(f,d);
  fclose(f);
  f = fopen('spektrum2_h.tab','w');
  fwrite(f,t);
  fclose(f);
