%This function estimates b_b and a on the basis of a given signal
%Detecttot(r), for rmax=3.1cm, and some detector geometry specifications 
%(deltar, Nshell). 
%This program uses the external function Energy_analysis.m. 
%Make sure this file is in the current dirrectory!
%For a 3.1 cm-radius detector with a radial resolution of 1cm, use
%lowres_Nshell=4 and lowres_deltar=(0.031-0.001)/(Nshell-1)=0.01.

function [bb, a] = Inversion_algorithm(Detectr, Detecttot, ...
                   lowres_deltar, lowres_Nshell)
%%%%%%%%%%%%%%%%%%%
%Geometry analysis%
%%%%%%%%%%%%%%%%%%%

lowres_Detecttot(1:lowres_Nshell)=0; %initialize low resolution
                                                 %Detecttot
deltar=0.001; %deltar is the 1mm-resolution used to
              %store the simulated backscattered
              %signal - see Monte Carlo simulation
              %code
lowres_Detectr(1:lowres_Nshell)=NaN;
lowres_Detectr(1)=0.001/2; %position associated to central 1mm disc
lowres_Detectr(2:lowres_Nshell) = ...
 (0.001+lowres_deltar/2):lowres_deltar:(0.001+(lowres_Nshell-1.5)*...
 lowres_deltar);              
for i=2:lowres_Nshell; %exclude first 1mm around the center from the 
                       %analysis
    for j=1:(lowres_deltar/deltar); %deltar is the 1mm-resolution used to
                                    %store the simulated backscattered
                                    %signal - see Monte Carlo simulation
                                    %code
        %Transform Monte Carlo model output (1mm-resolution) into simulated 
        %low-resolution instrument response needed for the IOP inversion: 
        lowres_Detecttot(i)=lowres_Detecttot(i)+...
            Detecttot(1+(i-2)*lowres_deltar/deltar+j);        
    end;    
end;    

Cumulative_Detecttot(1)=lowres_Detecttot(1); 
%lowres_Detecttot was already calculated excluding the first 1mm
for k=2:lowres_Nshell;
      Cumulative_Detecttot(k)=Cumulative_Detecttot(k-1)+...
                                        lowres_Detecttot(k);
end;    
%Exclude first value of Cumulative_Detecttot (and corresponding value of 
%lowres_Detectr) from the analysis, as the central 0 value should not 
%influence the fit:
fit_C(1:lowres_Nshell-1)=NaN;
fit_Detectr(1:lowres_Nshell-1)=NaN;
for i=1:(lowres_Nshell-1);
    fit_C(i)=Cumulative_Detecttot(i+1);
    fit_Detectr(i)=lowres_Detectr(i+1);
end;    
p2_1=polyfit(fit_Detectr,fit_C(:)',1);%use for linear fit
            alpha=p2_1(1);
            beta=p2_1(2);

%%%%%%%%%%%%%%%%%%%

[Detect_total]=Energy_analysis(Detectr, Detecttot);

bb=10^(1.048*log10(alpha)+0.3409);
bbp=bb-4.0883*10^(-4);
a=bb/(10^(-0.07407*(log10(Detect_total))^2+0.3525...
  *log10(Detect_total)+0.6555));
ap=a-0.43;

 disp([' alpha = ' num2str(alpha) 'm^{-1}']);
 disp([' Detect_total = ' num2str(Detect_total)]);
 disp([' bb = ' num2str(bb) 'm^{-1}']);
 disp([' bbp = ' num2str(bbp) 'm^{-1}']);
 disp([' a = ' num2str(a) 'm^{-1}']);
 disp([' a_p = ' num2str(ap) 'm^{-1}']);
      
return 

