%program solving Lab 4 homework

conc=[0 0.2404 0.4808 0.7212 1.2019 1.6827 2.4038 3.3654 5.0481];
count=[64.3706 178.0725 183.8979 186.5151 187.8622 186.8362 184.4440 178.7486 168.8259];

%Fitting the data, setting an over constrained linear matrix
%x_1(count_2-count_1)+x_2*(2*R*(C2-C1))=10log(C2/C1)
%where x_1 and x_2 are unknowns, C are concentrations, R the pathlength (5cm).
%1 and 2 denote two compared runs
%
% ignoring the first run with no beads (poor S/N) we have 7+6+5+..1=28 pairs

plot(conc(2:9),count(2:9));
xlabel('concentration [g/L]');
ylabel('ADV counts');

k=1
for i=3:9
    for j=2:i-1
    diff_count(k)=count(i)-count(j);
    diff_conc(k)=conc(i)-conc(j);
    conc_ratio(k)=conc(i)/conc(j);
    k=k+1;
    end
end
%setting the inversion problem. Since it is over constrained it will give
%us the best solution in a least-square sense.
% The vector x contain the answer we are looking for.

b=10*log10(conc_ratio)';

A=[diff_count; 2*0.05*diff_conc]';
[Q R]=qr(A);
x = R\(R'\(A'*b));
r = b - A*x;
e = R\(R'\(A'*r));
x = x + e;

%now solving the full problem:
%x1(3) x counts=K1-2 R x1(1)*C +10log10(x1(2) x C + K2)
%K1 = x1(3) * count(air), K2=10^(((count(blank)-count(air)) x x1(3))
%rewrite as: counts=60-2 R y1(1)*C+y1(2)10log10(y1(3) x C + 10^(0.43706/ y1(2)))

opts = optimset('fminsearch');      
opts = optimset(opts,'MaxIter',8000); 
opts = optimset(opts,'MaxFunEvals',4000);   % usually 100*number of params
opts = optimset(opts,'TolFun',1e-9);
%guess for paramters (use previous solutions plus a guess)
count_air=60;
count_blank=64.37;
x0=[x(2) 1/x(1) x(2)/4/pi];

%minimization routine
y1 = fminsearch(@least_square_backscattering,x0,opts,conc(1:9),count(1:9),count_air,count_blank)

%plot data and fit
R=0.05;
diff_count=count_blank-count_air;
plot(conc(1:9),count(1:9),'*b',conc(1:9),count_air-2*R*y1(1)*conc+y1(2)*10*log10(y1(3)*conc+10^(diff_count/10/y1(2))),'k')

%compare both solution methods: values displayed are multiplicative
%constant [dB/counts] and alpha_p* [db/m]

x

[1/y1(2), y1(1)/y1(2)]


