function [a, b, SE_a, SE_b, cov_squared, Struct_fit] = Calibrate3Days_Powerfit(Struct, Avg, SE)
% [a, b, SE_a, SE_b, cov_squared, Struct_fit] = PowerfitRegression(Struct, Avg, SE)
% Powefit (log10) regression between the measurement of each OPC (as x),
% the mean of 5 OPCs (as y) with standard deviation of 5 OPCs (as y_err).
% Obtain intercept a, slope b, error of intercept a, error of intercept b
% and covariance squared between a and b
%
% Scripts used: Linearfit.m
for i = 1:8  % 8 bins
    y = Avg(:,i); % mean of 5 OPCs, for i^{th} bin
    y_err = SE(:,i); 
    y_power = log10(y);
    y_err_power = sqrt((y_err./y/log(10)).^2); % error propagation
    
    for j = 1:5
        x = Struct.(['D',num2str(j)]).(['n',num2str(i)]); %linear measurements of each OPC, in #/m^3
        x_power = log10(x);
        [a(i,j), b(i,j), SE_a(i,j), SE_b(i,j), cov_squared(i,j)] = Linearfit(x_power, y_power, y_err_power);
        % use the function Linearfit.m in the same folder      
        
        y_powerfit = a(i,j) + x_power*b(i,j);
        y_err_powerfit = sqrt(SE_a(i,j)^2 + SE_b(i,j)^2*x_power.^2 + 2*cov_squared(i,j)*x_power);
        y_fit = 10.^y_powerfit;
        y_err_fit = sqrt((10.^y_powerfit*log(10).*y_err_powerfit).^2);
        
        Struct_fit.(['D',num2str(j)]).(['n',num2str(i)]).y_fit = y_fit;
        Struct_fit.(['D',num2str(j)]).(['n',num2str(i)]).y_err_fit = y_err_fit;
        Struct_fit.(['D',num2str(j)]).(['n',num2str(i)]).unit = '#/m^3 per 30 minutes';
    end   
    
end
display('%--- function Calibrate3Days_Powerfit() complete ----%');

end