%This algorithm uses the external procedures Geometry_analysis.m 
%and Energy_analysis.m. Make sure these files exist in the same 
%directory as the main code!! 
%The algorithm also loads several model output datasets. Make sure the 
%addresses for these output files are correct!

clear all;

%Define new low-resolution detector shell position vector:
lowres_Nshell = 4; %1 central 1mm-radius ring we will not account for in 
                   %the analysis (corresponds to the location of the light 
                   %source) + 3 1-cm thick detectors. 
lowres_deltar = (0.031-0.001)/(lowres_Nshell-1); %obtain thickness of 
                                                 %detector rings from the 
                                                 %maximum detector radius 
                                                 %of 3.1 cm, excluding the 
                                                 %first 1mm around the 
                                                 %center.
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);

%% Data initialization
%%
%%%%%%%%%%%%%%%%%%
%for Bp = 0.025: %
%%%%%%%%%%%%%%%%%%

load('exp46.mat');

exp46_Detecttot=Detecttot;
exp46_Detectnorm=Detectnorm;

Bp_025_Detecttot=Detecttot;
Bp_025_Detectnorm=Detectnorm;

Detect_total(1:Na,1:Nc)=NaN;
rhalf(1:Na,1:Nc)=-1; %dummy - we won't need this; we just kept it so 
                     %we didn't have to change the old structure of 
                     %Energy_analysis.m. Ideally, one day we would 
                     %properly clean up the codes...
for ia=1:Na;
    a=avec(ia);
    for ic=1:Nc;
        c=cvec(ia,ic);
        if(c>a);
            [Detect_total(ia,ic)]= ...
             Energy_analysis(Detectr,Bp_025_Detecttot(ia,ic,:));
        end;
    end;
end;

lowres_Detecttot(1:Na, 1:Nc, 1:lowres_Nshell)=0; %initialize low resolution
                                                 %Detecttot
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)+...
            Bp_025_Detecttot(:,:,1+(i-2)*lowres_deltar/deltar+j);        
    end;    
end;    
Geometry_analysis; %calculates the slope of a linear function 
                              %that is used to fit Cumulative_Detecttot=
                              %f(r), within a region close to the center 
                              %(3.1cm).
Bp_025_lowres_Detecttot=lowres_Detecttot;


Bp_025_alpha_vector = alpha1vector;
Bp_025_beta_vector = beta1vector;
bbp=bpvec*Bp;
Bp_025_bbp_vector(1:Na*Nc)=NaN;
Bp_025_Detect_total_vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        Bp_025_Detect_total_vector((i-1)*Nc+j)=Detect_total(i,j); 
        Bp_025_bbp_vector((i-1)*Nc+j)=bbp(i,j); 
    end;
end;


%define new colors for upcoming plots
o2 = [1 0.8 0];

%%
%%%%%%%%%%%%%%%%%
%for Bp = 0.02: %
%%%%%%%%%%%%%%%%%

load('exp26.mat');

exp26_Detecttot=Detecttot;
exp26_Detectnorm=Detectnorm;

load('exp28.mat');

exp28_Detecttot=Detecttot;
exp28_Detectnorm=Detectnorm;

load('exp29.mat');

exp29_Detecttot=Detecttot;
exp29_Detectnorm=Detectnorm;

load('exp36.mat');

exp36_Detecttot=Detecttot;
exp36_Detectnorm=Detectnorm;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%average out equivalent datasets%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Bp_02_Detecttot=(exp26_Detecttot+exp28_Detecttot+exp29_Detecttot+...
                exp36_Detecttot)/4;
Bp_02_Detectnorm=(exp26_Detectnorm+exp28_Detectnorm+exp29_Detectnorm...
                 +exp36_Detectnorm)/4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Detect_total(1:Na,1:Nc)=NaN;
rhalf(1:Na,1:Nc)=-1; %dummy 
for ia=1:Na;
    a=avec(ia);
    for ic=1:Nc;
        c=cvec(ia,ic);
        if(c>a);
            [Detect_total(ia,ic)]= ...
             Energy_analysis(Detectr,Bp_02_Detecttot(ia,ic,:));
        end;
    end;
end;

lowres_Detecttot(1:Na, 1:Nc, 1:lowres_Nshell)=0; %initialize low resolution
                                                 %Detecttot
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)+...
            Bp_02_Detecttot(:,:,1+(i-2)*lowres_deltar/deltar+j);        
    end;    
end;    
Geometry_analysis; %calculates the slope of a linear function 
                              %that is used to fit Cumulative_Detecttot=
                              %f(r), within a region close to the center 
                              %(3.1cm).
Bp_02_lowres_Detecttot=lowres_Detecttot;

Bp_02_alpha_vector = alpha1vector;
Bp_02_beta_vector = beta1vector;
bbp=bpvec*Bp;
Bp_02_bbp_vector(1:Na*Nc)=NaN;
Bp_02_Detect_total_vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        Bp_02_Detect_total_vector((i-1)*Nc+j)=Detect_total(i,j); 
        Bp_02_bbp_vector((i-1)*Nc+j)=bbp(i,j); 
    end;
end;

%%
%%%%%%%%%%%%%%%%%%
%for Bp = 0.015: %
%%%%%%%%%%%%%%%%%%

load('exp47.mat');

exp47_Detecttot=Detecttot;
exp47_Detectnorm=Detectnorm;

Bp_015_Detecttot=Detecttot;
Bp_015_Detectnorm=Detectnorm;


Detect_total(1:Na,1:Nc)=NaN;
rhalf(1:Na,1:Nc)=-1; %dummy
for ia=1:Na;
    a=avec(ia);
    for ic=1:Nc;
        c=cvec(ia,ic);
        if(c>a);
            [Detect_total(ia,ic)]= ...
             Energy_analysis(Detectr,Bp_015_Detecttot(ia,ic,:));
        end;
    end;
end;

lowres_Detecttot(1:Na, 1:Nc, 1:lowres_Nshell)=0; %initialize low resolution
                                                 %Detecttot
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)+...
            Bp_015_Detecttot(:,:,1+(i-2)*lowres_deltar/deltar+j);        
    end;    
end;    
Geometry_analysis; %calculates the slope of a linear function 
                              %that is used to fit Cumulative_Detecttot=
                              %f(r), within a region close to the center 
                              %(3.1cm).
Bp_015_lowres_Detecttot=lowres_Detecttot;

Bp_015_alpha_vector = alpha1vector;
Bp_015_beta_vector = beta1vector;
bbp=bpvec*Bp;
Bp_015_bbp_vector(1:Na*Nc)=NaN;
Bp_015_Detect_total_vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        Bp_015_Detect_total_vector((i-1)*Nc+j)=Detect_total(i,j); 
        Bp_015_bbp_vector((i-1)*Nc+j)=bbp(i,j); 
    end;
end;

%define new colors for upcoming plots
m3 = [0.6 0.2 0.8];

%%
%%%%%%%%%%%%%%%%%
%for Bp = 0.01: %
%%%%%%%%%%%%%%%%%

load('exp31.mat');

exp31_Detecttot=Detecttot;
exp31_Detectnorm=Detectnorm;


load('exp33.mat');

exp33_Detecttot=Detecttot;
exp33_Detectnorm=Detectnorm;


load('exp35.mat');

exp35_Detecttot=Detecttot;
exp35_Detectnorm=Detectnorm;


load('exp38.mat');

exp38_Detecttot=Detecttot;
exp38_Detectnorm=Detectnorm;


load('exp39.mat');

exp39_Detecttot=Detecttot;
exp39_Detectnorm=Detectnorm;


load('exp45.mat');

exp45_Detecttot=Detecttot;
exp45_Detectnorm=Detectnorm;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%average out equivalent datasets%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Bp_01_Detecttot=(exp31_Detecttot+exp33_Detecttot+exp35_Detecttot+...
                exp38_Detecttot+exp39_Detecttot+exp45_Detecttot)/6;
Bp_01_Detectnorm=(exp31_Detectnorm+exp33_Detectnorm+...
                 exp35_Detectnorm+exp38_Detectnorm+exp39_Detectnorm+...
                 exp45_Detectnorm)/6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Detect_total(1:Na,1:Nc)=NaN;
rhalf(1:Na,1:Nc)=-1; %dummy 
for ia=1:Na;
    a=avec(ia);
    for ic=1:Nc;
        c=cvec(ia,ic);
        if(c>a);
            [Detect_total(ia,ic)]= ...
             Energy_analysis(Detectr,Bp_01_Detecttot(ia,ic,:));
        end;
    end;
end;

lowres_Detecttot(1:Na, 1:Nc, 1:lowres_Nshell)=0; %initialize low resolution
                                                 %Detecttot
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)+...
            Bp_01_Detecttot(:,:,1+(i-2)*lowres_deltar/deltar+j);        
    end;    
end;    
Geometry_analysis; %calculates the slope of a linear function 
                              %that is used to fit Cumulative_Detecttot=
                              %f(r), within a region close to the center 
                              %(3.1cm).
Bp_01_lowres_Detecttot=lowres_Detecttot;

Bp_01_alpha_vector = alpha1vector;
Bp_01_beta_vector = beta1vector;
bbp=bpvec*Bp;
Bp_01_bbp_vector(1:Na*Nc)=NaN;
Bp_01_Detect_total_vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        Bp_01_Detect_total_vector((i-1)*Nc+j)=Detect_total(i,j); 
        Bp_01_bbp_vector((i-1)*Nc+j)=bbp(i,j); 
    end;
end;

%%
%%%%%%%%%%%%%%%%%%
%for Bp = 0.005: %
%%%%%%%%%%%%%%%%%%

load('exp30.mat');

exp30_Detecttot=Detecttot;
exp30_Detectnorm=Detectnorm;


load('exp32.mat');

exp32_Detecttot=Detecttot;
exp32_Detectnorm=Detectnorm;


load('exp34.mat');

exp34_Detecttot=Detecttot;
exp34_Detectnorm=Detectnorm;


load('exp40.mat');

exp40_Detecttot=Detecttot;
exp40_Detectnorm=Detectnorm;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%average out equivalent datasets%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Bp_005_Detecttot=(exp30_Detecttot+exp32_Detecttot+exp34_Detecttot+...
                 exp40_Detecttot)/4;
Bp_005_Detectnorm=(exp30_Detectnorm+exp32_Detectnorm+...
                  exp34_Detectnorm+exp40_Detectnorm)/4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Detect_total(1:Na,1:Nc)=NaN;
rhalf(1:Na,1:Nc)=-1; %dummy
for ia=1:Na;
    a=avec(ia);
    for ic=1:Nc;
        c=cvec(ia,ic);
        if(c>a);
            [Detect_total(ia,ic)]= ...
             Energy_analysis(Detectr,Bp_005_Detecttot(ia,ic,:));
        end;
    end;
end;

lowres_Detecttot(1:Na, 1:Nc, 1:lowres_Nshell)=0; %initialize low resolution
                                                 %Detecttot
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)+...
            Bp_005_Detecttot(:,:,1+(i-2)*lowres_deltar/deltar+j);        
    end;    
end;    
Geometry_analysis; %calculates the slope of a linear function 
                              %that is used to fit Cumulative_Detecttot=
                              %f(r), within a region close to the center 
                              %(3.1cm).
Bp_005_lowres_Detecttot=lowres_Detecttot;

Bp_005_alpha_vector = alpha1vector;
Bp_005_beta_vector = beta1vector;
bbp=bpvec*Bp;
Bp_005_bbp_vector(1:Na*Nc)=NaN;
Bp_005_Detect_total_vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        Bp_005_Detect_total_vector((i-1)*Nc+j)=Detect_total(i,j); 
        Bp_005_bbp_vector((i-1)*Nc+j)=bbp(i,j); 
    end;
end;


%For the graphical representations of the results, we need to convert 
%our analysis data to vector format:
avector(1:Na*Nc)=NaN;
cvector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        avector((i-1)*Nc+j)=avec(i); 
    end;
end;
for i=1:Na; 
    for j=1:Nc; 
        cvector((i-1)*Nc+j)=cvec(i,j); 
    end;

end;

%% Error analysis
%%%%%%%%%%%%%%%%
%generate plots%
%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Put all the data together%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

bbp_all(1:100)=Bp_025_bbp_vector(1:100);
alpha_all(1:100)=Bp_025_alpha_vector(1:100);
bbp_all(101:200)=Bp_02_bbp_vector(1:100);
alpha_all(101:200)=Bp_02_alpha_vector(1:100);
bbp_all(201:300)=Bp_015_bbp_vector(1:100);
alpha_all(201:300)=Bp_015_alpha_vector(1:100);
bbp_all(301:400)=Bp_01_bbp_vector(1:100);
alpha_all(301:400)=Bp_01_alpha_vector(1:100);
bbp_all(401:500)=Bp_005_bbp_vector(1:100);
alpha_all(401:500)=Bp_005_alpha_vector(1:100);

%%histogram of bbp values
%figure(1);
%hist(bbp_all,50);
%xlabel('b_{bp} [1/m]');
%title(['Histogram of b_{bp} values for the combined data sets of '...
%       'experiments Bp_{-}025, Bp_{-}02, Bp_{-}015, Bp_{-}01, '...
%        'Bp_{-}005']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Check relationship between bb_total and alpha%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Bp_025_bb_vector=Bp_025_bbp_vector+0.00040883;
                 %b_b,sw=4.0883*10^(-4)m^(-1) 
                 %- see Alina's notes/p.2/April 7, 2010
Bp_02_bb_vector=Bp_02_bbp_vector+0.00040883;
Bp_015_bb_vector=Bp_015_bbp_vector+0.00040883;
Bp_01_bb_vector=Bp_01_bbp_vector+0.00040883;
Bp_005_bb_vector=Bp_005_bbp_vector+0.00040883;
bb_all=bbp_all+0.00040883;

% We cannot use polyfit to fit bb vs. alpha because of the NaN values.
% Interactively, use cftool to obtain fits
%figure(2);
%plot(Bp_025_bb_vector,-(Bp_025_bb_vector-10.^(1.048*...
%     log10(Bp_025_alpha_vector)+0.3409))./Bp_025_bb_vector,'x',...
%     'Color',o2);
%hold on;
%plot(Bp_02_bb_vector,-(Bp_02_bb_vector-10.^(1.048*...
%     log10(Bp_02_alpha_vector)+0.3409))./Bp_02_bb_vector,'rx');
%plot(Bp_015_bb_vector,-(Bp_015_bb_vector-10.^(1.048*...
%     log10(Bp_015_alpha_vector)+0.3409))./Bp_015_bb_vector,'x',...
%     'Color',m3);
%plot(Bp_01_bb_vector,-(Bp_01_bb_vector-10.^(1.048*...
%     log10(Bp_01_alpha_vector)+0.3409))./Bp_01_bb_vector,'bx');
%plot(Bp_005_bb_vector,-(Bp_005_bb_vector-10.^(1.048*...
%     log10(Bp_005_alpha_vector)+0.3409))./Bp_005_bb_vector,'gx');
%xlabel('b_{b} [1/m]');
%ylabel('(10^{1.048*log10(alpha)+0.3409}-b_{b})/b_{b} ');
%legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
%       'EastOutside');
%hold off;

error_bb_all=(10.^(1.048*log10(alpha_all)+0.3409)-bb_all)./bb_all;
figure(3);
hist(100*error_bb_all,50);
xlabel('b_{b} log-log fit inversion relative errors [%]');
title(['Histogram for the combined data sets of experiments ' ...
       'Bp_{-}025, Bp_{-}02, Bp_{-}015, Bp_{-}01, Bp_{-}005']);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Get relative error ranges for 90% of the data%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sorted_error_bb_all=sort(abs(error_bb_all));
sorted_error_bb_all(isnan(sorted_error_bb_all))=[];
[aux,aux_size]=size(sorted_error_bb_all);
sample_error_bb(1:ceil(0.9*aux_size))=NaN;%this vector will contain 
                                           %the bb relative error 
                                           %values for the 90% of the 
                                           %data with the lowest errors
                                           %in bb inversion.
for i=1:ceil(0.9*aux_size);
    sample_error_bb(i)=sorted_error_bb_all(i);
end;

display(['90% of the b_b values are inverted with an accuracy ' ...
         'within ' sprintf('%g',100*max(abs(sample_error_bb))) '%']);
display(['100% of the b_b values are inverted with an accuracy ' ...
         'within ' sprintf('%g',100*max(abs(error_bb_all))) '%']);
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Get relative error ranges for bbp for 90% of the data%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
error_bbp_all=((10.^(1.048*log10(alpha_all)+0.3409)-...
              4.0883*10^(-4))-bbp_all)./bbp_all;
%figure(4);
%hist(100*error_bbp_all,50);
%xlabel('b_{bp} log-log fit inversion relative errors [%]');
%title(['Histogram for the combined data sets of experiments' ...
%    'Bp_{-}025, Bp_{-}02, Bp_{-}015, Bp_{-}01, Bp_{-}005']);

sorted_error_bbp_all=sort(abs(error_bbp_all));
sorted_error_bbp_all(isnan(sorted_error_bbp_all))=[];
[aux,aux_size]=size(sorted_error_bbp_all);
sample_error_bbp(1:ceil(0.9*aux_size))=NaN;%this vector will contain 
                                           %the bbp relative error 
                                           %values for the 90% of the 
                                           %data with the lowest errors
                                           %in bbp inversion.
for i=1:ceil(0.9*aux_size);
    sample_error_bbp(i)=sorted_error_bbp_all(i);
end;

display(['90% of the b_bp values are inverted with an accuracy ' ...
         'within ' sprintf('%g',100*max(abs(sample_error_bbp))) '%']);
display(['100% of the b_bp values are inverted with an accuracy ' ...
         'within ' sprintf('%g',100*max(abs(error_bbp_all))) '%']);

Bp_025_bb_fit=10.^(1.048*log10(Bp_025_alpha_vector)+0.3409);
Bp_02_bb_fit=10.^(1.048*log10(Bp_02_alpha_vector)+0.3409);
Bp_015_bb_fit=10.^(1.048*log10(Bp_015_alpha_vector)+0.3409);
Bp_01_bb_fit=10.^(1.048*log10(Bp_01_alpha_vector)+0.3409);
Bp_005_bb_fit=10.^(1.048*log10(Bp_005_alpha_vector)+0.3409);
inverted_bb_all=10.^(1.048*log10(alpha_all)+0.3409);
alongvector=repmat(avector,1,5);%vector of a values to match bb_all
Detect_total_all(1:100)=Bp_025_Detect_total_vector;
Detect_total_all(101:200)=Bp_02_Detect_total_vector;
Detect_total_all(201:300)=Bp_015_Detect_total_vector;
Detect_total_all(301:400)=Bp_01_Detect_total_vector;
Detect_total_all(401:500)=Bp_005_Detect_total_vector;

%%%%%%%%%%%%%
%a inversion%
%%%%%%%%%%%%%

%Log-log quadratic fit (obtained with cftool): 
%bb/a=10^{-0.07407*(log10(Detect_total))^2+0.3525* ...
%log10(Detect_total)+0.6555}
Bp_025_bb_a_fit=10.^(-0.07407*(log10(Bp_025_Detect_total_vector)) ...
                .^2+0.3525*log10(Bp_025_Detect_total_vector)+0.6555);
Bp_02_bb_a_fit=10.^(-0.07407*(log10(Bp_02_Detect_total_vector)) ...
               .^2+0.3525*log10(Bp_02_Detect_total_vector)+0.6555);
Bp_015_bb_a_fit=10.^(-0.07407*(log10(Bp_015_Detect_total_vector)) ...
                .^2+0.3525*log10(Bp_015_Detect_total_vector)+0.6555);
Bp_01_bb_a_fit=10.^(-0.07407*(log10(Bp_01_Detect_total_vector)) ...
               .^2+0.3525*log10(Bp_01_Detect_total_vector)+0.6555);
Bp_005_bb_a_fit=10.^(-0.07407*(log10(Bp_005_Detect_total_vector)) ...
                .^2+0.3525*log10(Bp_005_Detect_total_vector)+0.6555);

%figure(5);
%semilogx(avector,-(avector-(Bp_025_bb_fit./Bp_025_bb_a_fit)) ...
%./avector,'x','Color',o2);
%hold on;
%semilogx(avector,-(avector-(Bp_02_bb_fit./Bp_02_bb_a_fit)) ...
%./avector,'rx');
%semilogx(avector,-(avector-(Bp_015_bb_fit./Bp_015_bb_a_fit)) ...
%./avector,'x','Color',m3);
%semilogx(avector,-(avector-(Bp_01_bb_fit./Bp_01_bb_a_fit)) ...
%./avector,'bx');
%semilogx(avector,-(avector-(Bp_005_bb_fit./Bp_005_bb_a_fit)) ...
%./avector,'gx');
%title('b_b-based a inversion');
%xlabel('a [1/m]');
%ylabel('(a_{inverted}-a_{empirical})/a_{empirical} ');
%legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location', ...
%       'EastOutside');
%hold off;

%histogram of relative errors in a inversion using the bb-based fit

error_a_all(1:100)=-(avector-(Bp_025_bb_fit./Bp_025_bb_a_fit)) ...
                   ./avector;
error_a_all(101:200)=-(avector-(Bp_02_bb_fit./Bp_02_bb_a_fit)) ...
                     ./avector;
error_a_all(201:300)=-(avector-(Bp_015_bb_fit./Bp_015_bb_a_fit)) ...
                     ./avector;
error_a_all(301:400)=-(avector-(Bp_01_bb_fit./Bp_01_bb_a_fit)) ...
                     ./avector;
error_a_all(401:500)=-(avector-(Bp_005_bb_fit./Bp_005_bb_a_fit)) ...
                     ./avector;

figure(6);
hist(100*error_a_all,50);
xlabel('a inversion relative errors [%]');
title(['Histogram for the combined data sets of experiments ' ...
       'Bp_{-}025, Bp_{-}02, Bp_{-}015, Bp_{-}01, Bp_{-}005']);

sorted_error_a_all=sort(abs(error_a_all));
sorted_error_a_all(isnan(sorted_error_a_all))=[];
[aux,aux_size]=size(sorted_error_a_all);
sample_error_a(1:ceil(0.9*aux_size))=NaN;%this vector will contain 
                                         %the a relative error values 
                                         %for the 90% of the data with 
                                         %the lowest errors in bbp 
                                         %inversion.
for i=1:ceil(0.9*aux_size);
    sample_error_a(i)=sorted_error_a_all(i);
end;

display(['90% of the a values are inverted with and accuracy ' ...
         'within ' sprintf('%g',100*max(abs(sample_error_a))) '%.']);
display(['100% of the a values are inverted with and accuracy ' ...
         'within ' sprintf('%g',100*max(abs(error_a_all))) '%.']);
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate errors for ap inversion%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%histogram of relative errors in a inversion using the bb-based fit

error_ap_all(1:100)=-(avector-(Bp_025_bb_fit./Bp_025_bb_a_fit)) ...
                   ./(avector-asw);
error_ap_all(101:200)=-(avector-(Bp_02_bb_fit./Bp_02_bb_a_fit)) ...
                     ./(avector-asw);
error_ap_all(201:300)=-(avector-(Bp_015_bb_fit./Bp_015_bb_a_fit)) ...
                     ./(avector-asw);
error_ap_all(301:400)=-(avector-(Bp_01_bb_fit./Bp_01_bb_a_fit)) ...
                     ./(avector-asw);
error_ap_all(401:500)=-(avector-(Bp_005_bb_fit./Bp_005_bb_a_fit)) ...
                     ./(avector-asw);

%figure(7);
%hist(100*error_ap_all,50);
%xlabel('a_p inversion relative errors [%]');
%title(['Histogram for the combined data sets of experiments ' ...
%       'Bp_{-}025, Bp_{-}02, Bp_{-}015, Bp_{-}01, Bp_{-}005']);

sorted_error_ap_all=sort(abs(error_ap_all));
sorted_error_ap_all(isnan(sorted_error_ap_all))=[];
[aux,aux_size]=size(sorted_error_ap_all);
sample_error_ap(1:ceil(0.9*aux_size))=NaN;%this vector will contain 
                                          %the ap relative error values 
                                          %for the 90% of the data with 
                                          %the lowest errors in bbp 
                                          %inversion.
for i=1:ceil(0.9*aux_size);
    sample_error_ap(i)=sorted_error_ap_all(i);
end;

display(['90% of the a_p values are inverted with and accuracy ' ...
         'within ' sprintf('%g',100*max(abs(sample_error_ap))) '%.']);
display(['100% of the a_p values are inverted with and accuracy ' ...
         'within ' sprintf('%g',100*max(abs(error_ap_all))) '%.']);
     
%% Explore distribution of inversion errors
Bp_025_a_errors=-(avector-(Bp_025_bb_fit./Bp_025_bb_a_fit))./avector;
Bp_02_a_errors=-(avector-(Bp_02_bb_fit./Bp_02_bb_a_fit))./avector;
Bp_015_a_errors=-(avector-(Bp_015_bb_fit./Bp_015_bb_a_fit))./avector;
Bp_01_a_errors=-(avector-(Bp_01_bb_fit./Bp_01_bb_a_fit))./avector;
Bp_005_a_errors=-(avector-(Bp_005_bb_fit./Bp_005_bb_a_fit))./avector;

% Relative errors in a inversion vs. bb/a value
figure(9);
semilogx(Bp_025_bb_vector./avector,100*Bp_025_a_errors,'x','Color',o2);
hold on;
semilogx(Bp_02_bb_vector./avector,100*Bp_02_a_errors,'rx');
semilogx(Bp_015_bb_vector./avector,100*Bp_015_a_errors,'x','Color',m3);
semilogx(Bp_01_bb_vector./avector,100*Bp_01_a_errors,'bx');
semilogx(Bp_005_bb_vector./avector,100*Bp_005_a_errors,'gx');
xlabel('b_b/a');
ylabel('a inversion relative errors [%] ');
legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
       'EastOutside');
hold off;

% Relative errors in a inversion vs. single scattering albedo value
figure(10);
plot((cvector-avector)./cvector,100*Bp_025_a_errors,'x','Color',o2);
hold on;
plot((cvector-avector)./cvector,100*Bp_02_a_errors,'rx');
plot((cvector-avector)./cvector,100*Bp_015_a_errors,'x','Color',m3);
plot((cvector-avector)./cvector,100*Bp_01_a_errors,'bx');
plot((cvector-avector)./cvector,100*Bp_005_a_errors,'gx');
xlabel('b/c');
ylabel('a inversion relative errors [%] ');
legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
       'EastOutside');
hold off;

%%Relative errors in a inversion vs. particle single scattering albedo value
%Bp_025_omegap=Bp_025_bbp_vector./(Bp_025_bbp_vector+0.025*(avector-asw));
%Bp_02_omegap=Bp_02_bbp_vector./(Bp_02_bbp_vector+0.02*(avector-asw));
%Bp_015_omegap=Bp_015_bbp_vector./(Bp_015_bbp_vector+0.015*(avector-asw));
%Bp_01_omegap=Bp_01_bbp_vector./(Bp_01_bbp_vector+0.01*(avector-asw));
%Bp_005_omegap=Bp_005_bbp_vector./(Bp_005_bbp_vector+0.005*(avector-asw));
%figure(13);
%plot(Bp_025_omegap,100*Bp_025_a_errors,'x','Color',o2);
%hold on;
%plot(Bp_02_omegap,100*Bp_02_a_errors,'rx');
%plot(Bp_015_omegap,100*Bp_015_a_errors,'x','Color',m3);
%plot(Bp_01_omegap,100*Bp_01_a_errors,'bx');
%plot(Bp_005_omegap,100*Bp_005_a_errors,'gx');
%xlabel('b_p/c_p');
%ylabel('a inversion relative errors [%] ');
%legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
%       'EastOutside');
%hold off;


Bp_025_bb_errors=-(Bp_025_bb_vector-Bp_025_bb_fit)./Bp_025_bb_vector;
Bp_02_bb_errors=-(Bp_02_bb_vector-Bp_02_bb_fit)./Bp_02_bb_vector;
Bp_015_bb_errors=-(Bp_015_bb_vector-Bp_015_bb_fit)./Bp_015_bb_vector;
Bp_01_bb_errors=-(Bp_01_bb_vector-Bp_01_bb_fit)./Bp_01_bb_vector;
Bp_005_bb_errors=-(Bp_005_bb_vector-Bp_005_bb_fit)./Bp_005_bb_vector;

% Relative errors in b_b inversion vs. bb/a value
figure(11);
semilogx(Bp_025_bb_vector./avector,100*Bp_025_bb_errors,'x','Color',o2);
hold on;
semilogx(Bp_02_bb_vector./avector,100*Bp_02_bb_errors,'rx');
semilogx(Bp_015_bb_vector./avector,100*Bp_015_bb_errors,'x','Color',m3);
semilogx(Bp_01_bb_vector./avector,100*Bp_01_bb_errors,'bx');
semilogx(Bp_005_bb_vector./avector,100*Bp_005_bb_errors,'gx');
xlabel('b_b/a');
ylabel('b_b inversion relative errors [%] ');
legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
       'EastOutside');
hold off;

% Relative errors in b_b inversion vs. single scattering albedo value
figure(12);
plot((cvector-avector)./cvector,100*Bp_025_bb_errors,'x','Color',o2);
hold on;
plot((cvector-avector)./cvector,100*Bp_02_bb_errors,'rx');
plot((cvector-avector)./cvector,100*Bp_015_bb_errors,'x','Color',m3);
plot((cvector-avector)./cvector,100*Bp_01_bb_errors,'bx');
plot((cvector-avector)./cvector,100*Bp_005_bb_errors,'gx');
xlabel('b/c');
ylabel('b_b inversion relative errors [%] ');
legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
       'EastOutside');
hold off;     

%figure(14);
%plot(Bp_025_omegap,100*Bp_025_bb_errors,'x','Color',o2);
%hold on;
%plot(Bp_02_omegap,100*Bp_02_bb_errors,'rx');
%plot(Bp_015_omegap,100*Bp_015_bb_errors,'x','Color',m3);
%plot(Bp_01_omegap,100*Bp_01_bb_errors,'bx');
%plot(Bp_005_omegap,100*Bp_005_bb_errors,'gx');
%xlabel('b_p/c_p');
%ylabel('b_b inversion relative errors [%] ');
%legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
%       'EastOutside');
%hold off;

%figure(15);
%plot(Bp_025_omegap,-100*(avector-(Bp_025_bb_fit./Bp_025_bb_a_fit)) ...
%                   ./(avector-asw),'x','Color',o2);
%hold on;
%plot(Bp_02_omegap,-100*(avector-(Bp_02_bb_fit./Bp_02_bb_a_fit)) ...
%                   ./(avector-asw),'rx');
%plot(Bp_015_omegap,-100*(avector-(Bp_015_bb_fit./Bp_015_bb_a_fit)) ...
%                   ./(avector-asw),'x','Color',m3);
%plot(Bp_01_omegap,-100*(avector-(Bp_01_bb_fit./Bp_01_bb_a_fit)) ...
%                   ./(avector-asw),'bx'); 
%plot(Bp_005_omegap,-100*(avector-(Bp_005_bb_fit./Bp_005_bb_a_fit)) ...
%                   ./(avector-asw),'gx');       
%xlabel('b_p/c_p');
%ylabel('a_p inversion relative errors [%] ');
%legend('Bp=2.5%','Bp=2%','Bp=1.5%','Bp=1%','Bp=0.5%','Location',...
%       'EastOutside');
%hold off;
               