clear all;
d=load('alu.dat');
r=d(:,1);
n=d(:,2)*10;
d=r;
dd=0.01;

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

%standartabweichung
dn=sqrt(n);


%Fit
%   f= y0 + A exp(-µ1 d) + B exp(-µ2 d)
%   f= y(1) + y(2) exp(-y(3) d) + y(4) exp(-y(5) d)
     a=[0,500,16,500,2];
 	eps = 0.0001;     %max relativer fehler
	
	ndata = numel(n);
	ma = numel(a);
	[a,sd_a,var,covar,c] = mrqmin(d,n,dn,ndata,ma,a,eps,'funcs_alu');
     disp('f= a(1) + a(2) exp(-a(3) d) + a(4) exp(-a(5) d)');
     disp(a);
     disp(sd_a);
    
        dd = linspace(0,max(d)*1.1);
	f= a(1) + a(2).* exp(-a(3) .* dd) + a(4).* exp(-a(5) .* dd);
  	
	chi2=var(1)*(numel(d)-ma);
    	disp('chi^2');
	disp(chi2);
    	disp('Varianz');
    	disp(var);
    	sigma_v=sqrt(2/(numel(d)-4));
    	disp('Standartabweichung der Varianz');
    	disp(sigma_v);

	rho=2.7;
	E01 = (a(3)/rho)^(1.39);
	E02 = (a(5)/rho)^(1.39);
	dE01 = 1.39 * (a(3)/rho)^(0.39)*sd_a(3);
	dE02 = 1.39 * (a(5)/rho)^(0.39)*sd_a(5);

	disp('EO1, dEO1');
	disp([E01,dE01]);	
	disp('EO2, dEO2');
	disp([E02,dE02]);	


hold on;
	p3=plot(dd,f);
	p1=errorbar(d,n,dn,'.');
	ax1=gca;
	set(ax1,'YLim',[0,max(f)*1.1]);
	set(ax1,'XLim',[-1,max(dd)*1.1]);
	set(ax1,'YScale','log');
	
	le1=legend([p1(2),p3],'n','Doppelexponentialfit',1);
	xl1=xlabel('Aluminium- Dicke d/[mm]');
	yl1=ylabel('Zählwert n');

	grid on;
hold off;

print -deps2 -tiff -r300 absorption


t=['$d/[mm]$ & $n$ & $dn$ \\'];
d=latex([d,n,dn],'%3.1f','nomath');
f = fopen('absorption.tab','w');
fwrite(f,d);
fclose(f);
f = fopen('absorption_h.tab','w');
fwrite(f,t);
fclose(f);


