clc
clear all 
close all

%%%INSERT THE DATA%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Insert here the current steps, the standard deviations, and the
%%concentrations
concentration=[];%!!!!!!!!!!!!!!!!!!!!!! in mM!!!!!
mean_avg=[];
dev_std=[];

% concentration=[];%!!!!!!!!!!!!!!!!!!!!!! in mM!!!!!
% mean_avg=[];
% dev_std=[];
% 

%blank measurements
b1=32.5016;
b2=32.3759;
b3=32.3183;
%STANDARD DEVIATION BALNK MEASUREMENT
std_blank=[0.0194,0.1378,0.0269];


num_concentrazioni=7; 
diameter_mm=1;
radius_mm=diameter_mm/2;
areael=pi.*radius_mm^2; %area electrode mm^2

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

%Evalutation of the mean and of standard deviation of the current for each
%concentration
% mean_avg=zeros(1,num_concentrazioni);
% dev_std=zeros(1,num_concentrazioni);
std_blank=zeros(1,3);

%1.
%Fitting and linear regression analysis (weighted with the standard
%deviations)

s = fitoptions('Weights', dev_std');
f = fittype('poly1');
[c2,gof2] = fit(concentration',mean_avg',f,s);

% c2 = 
% 
%      Linear model Poly1:
%      c2(x) = p1*x + p2
%      Coefficients (with 95% confidence bounds):
%        p1 =      0.4857  (0.4299, 0.5415)
%        p2 =       1.118  (-5.97, 8.206)

% 
% gof2 = 
% 
%            sse: 377.2462
%        rsquare: 0.9932
%            dfe: 4
%     adjrsquare: 0.9915
%           rmse: 9.7114

values=coeffvalues(c2)%on the main Matlab window I see the value for the intercept and the slope
ci = confint(c2,0.95); %it includes the upper and lower extremes of the confidence interval for slope and intercept
p12=predint(c2,concentration,0.95,'functional','on');%evaluation of the confidence intervals for each concentration
gof2

%Errors for slope and intercept
error_sensitivity=values(1,1)-ci(1,1)
error_intercept=values(1,2)-ci(1,2);

%Sensitivity normalization and relative error 
normalized_sensitivity=values(1,1)/areael
normalized_error_sensitivity=error_sensitivity/areael
relative_error_sensitivity=normalized_error_sensitivity/normalized_sensitivity;

%Residuals
yfit = c2(concentration)
yresid = mean_avg'-yfit;
SSresid = sum(yresid.^2);
SStotal = (length(c2(concentration)-1) * var(c2(concentration)));
rsq = 1 - SSresid/SStotal

% residual_std=sqrt(sum(yresid.^2)/(num_concentrazioni-2))

%Evaluation of the LOD
%LOD
mean_std_blank=mean(std_blank);
LOD=3*mean_std_blank/values(1,1)%here I needed the not-normalized sensitivity
%LOD2=3*residual_std/values(1,1)

% %Error on LOD
% std_std_blank=std(std_blank);
% relative_error_LOD=std_std_blank/mean_std_blank+relative_error_sensitivity;
% error_LOD=LOD*relative_error_LOD
%error_LOD2=LOD2*relative_error_LOD

%Final table with main data
TABLE_FINAL=[values,error_sensitivity,normalized_sensitivity,normalized_error_sensitivity,LOD];

%ANOVA table
lm = fitlm(concentration,mean_avg,'linear')
anova(lm,'summary')

%%%PLOTS%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calibration lines with error bars
figure()
errorbar(concentration,mean_avg,dev_std,'xr');
hold on
h1=plot(c2,'r');
hold on

title('Detection of Glucose 3-21mM');%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
xlabel('concentration [mM]');
ylabel('current [nA]');

%Calibration line with confidence interval
figure()
errorbar(concentration,mean_avg,dev_std,'xr');
title('Detection of Glucose 3-21mM');%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
hold on
plot(c2);
hold on
plot(concentration,p12,'m--');
xlabel('concentration [mM]');
ylabel('current [nA]');

%Residual Plot
figure()
plot(concentration,yresid,'--','Marker','o','MarkerFaceColor','red')
xlabel('concentration [mM]');
ylabel('Residuals');