clear all, close all
%
x=[1 3 13];
%y=[0.8 0.6 0.2]; % Hg
y=[79.7 63.2 18.2]; %Pb
%y=[208 145 109]; %Zn
%y=[2.8 0.6 0.1]; %Ag
%y=[0.5 0.3 0.2]; %Cd
%y=[175 161.4 136.4]; %V
%y=[40 27.6 14.1]; %Cu
%
figure(1), hold on

plot(0,0,'b-','linewidth',10)

rectangle('Position',[0 0 2 79.7],'facecolor','b')
rectangle('Position',[2 0 2 63.2],'facecolor','b')
rectangle('Position',[12 0 2 18.2],'facecolor','b')

x = x(:);
y = y(:);
%
st_ = [0.91170928542807805 -0.12118783131983744 ];
ft_ = fittype('exp1');
ok_ = isfinite(x) & isfinite(y);
%
% Fit this model using new data
cf = fit(x(ok_),y(ok_),ft_,'Startpoint',st_);
%
x1=linspace(0,30,200);
a=cf.a
b=cf.b
y1=a*exp(b*x1);
%
%plot(x1,y1,'r.-',x,y,'bo')
plot(x1,y1,'r-','LineWidth',2)
set(gca,'FontSize',14,'FontName','Arial', ...
        'XScale','lin','YScale','lin','box','on')
xlabel('Profondità [cm]')
ylabel('Concentrazione Pb [mg/kg]')
%
%
for ii=2:length(x1);
z1(ii)=trapz(x1(1:ii),y1(1:ii))/x1(ii);
end
%
zz1=a*(-1+exp(b*x1))/b./x1;
%
figure(1), hold on
%plot(x1,z1,'b.-',x1,zz1,'ro')
plot(x1,zz1,'k--','LineWidth',2)
set(gca,'FontSize',14,'FontName','Arial', ...
        'XScale','lin','YScale','lin','box','on')
xlabel('Profondità [cm]')
%ylabel('concentrazione media Pb [mg/kg]')

%
aa=sym('aa','real');
bb=sym('bb','real');
xx=sym('xx','real');
xxint=sym('xxint','real');
%
int(aa*exp(bb*xx),xx,0,xxint)
%
zz1(end)

legend('Misure Cibic et al.','Interpolazione','Valore medio')
legend boxoff