% For a detailed explanation of what the code does, check Alina's notes
% from page 4/Feb.10, 2010. 
% In short, for a given experiment that used code
% Monte_Carlo_RT_Model.m, this code 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 
% around the center, excluding the first 1mm). 

%load('exp_name.mat'); % you can load the data in the command window 
                  %too, or in a code that calls this as a subroutine


%define the matrices corresponding to the coefficients of the fit 
%function for Cumulative Detecttot vs. r, for each (a,c). After some 
%preliminary tests, we have decided to fit first-order polynomials, the
%coefficients of which will be alpha and beta.
alpha1(1:Na,1:Nc)=NaN;
beta1(1:Na,1:Nc)=NaN;

for i=1:Na;
    for j=1:Nc;
        Cumulative_Detecttot(i,j,1)=lowres_Detecttot(i,j,1); 
        %lowres_Detecttot was already calculated excluding the first 1mm
        for k=2:lowres_Nshell;
            Cumulative_Detecttot(i,j,k)=Cumulative_Detecttot(i,j,k-1)+...
                                        lowres_Detecttot(i,j,k);
        end;    
    end;
end;    

pC=permute(Cumulative_Detecttot,[3 1 2]);

%Exclude first value of pC (and corresponding value of lowres_Detectr) from
%the analysis, as the central 0 value should not influence the fit:
fit_pC(1:lowres_Nshell-1,1:Na,1:Nc)=NaN;
fit_Detectr(1:lowres_Nshell-1)=NaN;
for i=1:(lowres_Nshell-1);
    fit_pC(i,:,:)=pC(i+1,:,:);
    fit_Detectr(i)=lowres_Detectr(i+1);
end;    
    
for i=1:Na;
    for j=1:Nc;
        if (((i==1)&&(j==1))||(i~=1));
            p2_1=polyfit(fit_Detectr,fit_pC(:,i,j)',1);%use for linear fit
            alpha1(i,j)=p2_1(1);
            beta1(i,j)=p2_1(2);
        end;
    end;
end;    
            

%Create other vectors necessary in Data_analysis.m
cvector(1:Na*Nc)=0;
for i=1:Na; 
    for j=1:Nc; 
        cvector((i-1)*Nc+j)=cvec(i,j); 
    end;
end;
avector(1:Na*Nc)=0;
for i=1:Na; 
    for j=1:Nc; 
        avector((i-1)*Nc+j)=avec(i); 
    end;
end;
alpha1vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        alpha1vector((i-1)*Nc+j)=alpha1(i,j); 
    end;
end;
beta1vector(1:Na*Nc)=NaN;
for i=1:Na; 
    for j=1:Nc; 
        beta1vector((i-1)*Nc+j)=beta1(i,j); 
    end;
end;

return;