clc, clear, clf, close all; %clears cached data
binwidth=20;%number of points for each bin that will be fit
stepsize=10;% where the next range of points starts.if binwidth>stepsize, bins will overlap. if binwidth=stepsize, bins will not overlap.  Stepsize should not be greater than binwidth. Example: iff the binwidth is 15 and the step size is 10, the first bin is the first 15 points.  The second bin has 15 points, but the first 5 points are the last 5 points in the previous bin, and the next 10 points are "new" points in the fitting/binning
tol=1e-3; % sets the tolerances for how many points are included in the fit.  Raise to include more points, lower to decrease number of points
vidspeed=7452.68; % speed of the video, number of frame rates
data=[0	1
1.3418E-4	0.98579
2.6836E-4	0.97241
4.0254E-4	0.95903
5.3672E-4	0.94482
6.709E-4	0.9306
8.0508E-4	0.91639
9.39259E-4	0.90217
0.00107	0.88712
0.00121	0.87291
0.00134	0.85619
0.00148	0.84114
0.00161	0.82525
0.00174	0.80936
0.00188	0.79348
0.00201	0.77676
0.00215	0.76087
0.00228	0.74415
0.00242	0.72826
0.00255	0.7107
0.00268	0.69398
0.00282	0.67642
0.00295	0.65803
0.00309	0.64047
0.00322	0.62124
0.00335	0.60284
0.00349	0.58278
0.00362	0.56271
0.00376	0.54097
0.00389	0.5209
0.00403	0.49916
0.00416	0.47742
0.00429	0.45318
0.00443	0.4306
0.00456	0.40552
0.0047	0.38043
0.00483	0.35619
0.00496	0.3286
0.0051	0.30268
0.00523	0.27676
0.00537	0.25
0.0055	0.22575
0.00564	0.20234
0.00577	0.18562
0.0059	0.17391
0.00604	0.16806
0.00617	0.16388
0.00631	0.1597
0.00644	0.15803
0.00657	0.15635
0.00671	0.15385
0.00684	0.15217
0.00698	0.1505
0.00711	0.14799
0.00725	0.14632
0.00738	0.14465
0.00751	0.14214
0.00765	0.1388
0.00778	0.13629
0.00792	0.13462
0.00805	0.13211
0.00818	0.13043
0.00832	0.1296
0.00845	0.12709
0.00859	0.12542
0.00872	0.12458
0.00886	0.12207
0.00899	0.1204
0.00912	0.11957
0.00926	0.11874
0.00939	0.118
0.00953	0.11576
0.00966	0.11427
0.0098	0.11276
0.00993	0.11127
0.01006	0.11053
0.0102	0.10903
0.01033	0.10903
0.01047	0.10754
0.0106	0.10605
0.01073	0.10456
0.01087	0.10381
0.011	0.10306
0.01114	0.10231
0.01127	0.10007
0.01141	0.09932
0.01154	0.09858
0.01167	0.09783
0.01181	0.09709
0.01194	0.09559
0.01208	0.0941
0.01221	0.09336
0.01234	0.09261
0.01248	0.09261
0.01261	0.09111
0.01275	0.08961
0.01288	0.08887
0.01302	0.08812
0.01315	0.08812
0.01328	0.08738
0.01342	0.08514
0.01355	0.08588
0.01369	0.08439
0.01382	0.08365
0.01395	0.08215
0.01409	0.08215
0.01422	0.08065
0.01436	0.0799
0.01449	0.07916
0.01463	0.07841
0.01476	0.07767
0.01489	0.07617
0.01503	0.07468
0.01516	0.07468
0.0153	0.07468
0.01543	0.07394
0.01556	0.07394
0.0157	0.07319
0.01583	0.07319
0.01597	0.07244
0.0161	0.07095
0.01624	0.07095
0.01637	0.07095
0.0165	0.06945
0.01664	0.0687
0.01677	0.06796
0.01691	0.06721
0.01704	0.06646
0.01718	0.06646
0.01731	0.06572
0.01744	0.06497
0.01758	0.06497
0.01771	0.06423
0.01785	0.06423
0.01798	0.06273
0.01811	0.06348
0.01825	0.06124
0.01838	0.06199
0.01852	0.06124
0.01865	0.0605
0.01879	0.05975
0.01892	0.05975
0.01905	0.05899
0.01919	0.05825
0.01932	0.05825
0.01946	0.0575
0.01959	0.05601
0.01972	0.05601
0.01986	0.05601
0.01999	0.05526
0.02013	0.05452
0.02026	0.05377
0.0204	0.05377
0.02053	0.05228
0.02066	0.05228
0.0208	0.05153
0.02093	0.05153
0.02107	0.05079
0.0212	0.05004
0.02133	0.05004
0.02147	0.05004
0.0216	0.04929
0.02174	0.04779
0.02187	0.04705
0.02201	0.0463
0.02214	0.04705
0.02227	0.0463
0.02241	0.04555
0.02254	0.04555
0.02268	0.04406
0.02281	0.04406
0.02294	0.04406
0.02308	0.04331
0.02321	0.04331
0.02335	0.04182
0.02348	0.04182
0.02362	0.04182
0.02375	0.04108
0.02388	0.04033
0.02402	0.04108
0.02415	0.04033
0.02429	0.04033
0.02442	0.03884
0.02455	0.03884
0.02469	0.03809
0.02482	0.03659
0.02496	0.03659
0.02509	0.03734
0.02523	0.03584
0.02536	0.03584
0.02549	0.03584
0.02563	0.0351
0.02576	0.0351
0.0259	0.0351
0.02603	0.0351
0.02617	0.03435
0.0263	0.03435
0.02643	0.03361
0.02657	0.03361
0.0267	0.03361
0.02684	0.03286
0.02697	0.03286
0.0271	0.03211
0.02724	0.03211
0.02737	0.03211
0.02751	0.03137
0.02764	0.02987
0.02778	0.03062
0.02791	0.03062
0.02804	0.02913
0.02818	0.02913
0.02831	0.02913
0.02845	0.02838
0.02858	0.02838
0.02871	0.02838
0.02885	0.02764
0.02898	0.02764
0.02912	0.02764
0.02925	0.02688
0.02939	0.02613
0.02952	0.02613
0.02965	0.02613
0.02979	0.02539
0.02992	0.02539
0.03006	0.02613
0.03019	0.02464
0.03032	0.02464
0.03046	0.02464
0.03059	0.02464
0.03073	0.02464
0.03086	0.0239
0.031	0.0224
0.03113	0.0224
0.03126	0.02315
0.0314	0.02315
0.03153	0.0224
0.03167	0.02166
0.0318	0.02017
0.03193	0.02091
0.03207	0.02017
0.0322	0.01942
0.03234	0.01867
0.03247	0.01867
0.03261	0.01793
0.03274	0.01718
0.03287	0.01793
0.03301	0.01718
0.03314	0.01718
0.03328	0.01643
0.03341	0.01718
0.03354	0.01643
0.03368	0.01643
0.03381	0.01568
0.03395	0.01568
0.03408	0.01568
0.03422	0.01493
0.03435	0.01568
0.03448	0.01493
0.03462	0.01568
0.03475	0.01493
0.03489	0.01493
0.03502	0.01493
0.03516	0.01493
0.03529	0.01493
0.03542	0.01419
0.03556	0.01419
0.03569	0.01419
0.03583	0.01344
0.03596	0.01344
0.03609	0.01195
0.03623	0.01269
0.03636	0.01195
0.0365	0.0112
0.03663	0.0112
0.03677	0.01046
0.0369	0.00971
0.03703	0.00971
0.03717	0.00822
0.0373	0.00673
0.03744	0.00822
0.03757	0.00822
0.0377	0.00747
0.03784	0.00298
0.03797	0.00598
0.03811	0.00598
0.03824	0.00673
0.03838	0.00598
0.03851	0.00598
0.03864	0.00598
0.03878	0.00598
0.03891	0.00149]; %input radius by time data in a matrix with time data in column 1 and normalized radius data in column 2.  This can be directly copied from the Startpoint_Symmetry Code output, but must have an EC region to analyze.  An example dataset using 0.5 wt.% PEO is shown for reference

%%
%Normalizing data in case radius data is entere without normalizing
  xtot = data(:,1); % x=time
  ytot = data(:,2); %y=normalized radius
  maxytot = max(ytot); %finds the max radius value (should be first radius value at t=0)
  ytot = ytot/max(ytot);%normalizes so it starts at D/Do=1. If already normalized, will just divide 1/1
  endtime=xtot(end); %endtime is the total time taken for thinning 
  
xg=xtot; %matrix of time values for further processing, xtot is kept unprocessed and graphed as is
yg=ytot;%matrix of radius values for further processing, ytot is kept unprocessed and graphed as is

%Linear fit for each k bin, find the slope of logR vs. t data.  In the EC region, the slope should be constant (horizontal)
 
for k=1:(floor((length(xg)-(binwidth-stepsize))/stepsize)) %for loop for the kth bin where k goes from 1 to the total number of bins, subtracting binwidth-stepsize accounts for bin overlapping if stepsize<binwidth;
   xw=xg((1+stepsize*k-stepsize):(stepsize*k-stepsize+binwidth)); %determines the start and end time of bin k
   yw=yg((1+stepsize*k-stepsize):(stepsize*k-stepsize+binwidth));%determines the start and end radius of bin k
   linfit=fitlm(xw,log(yw)); % a linear fit over the log of radius and time value
   R_squared(k)=linfit.Rsquared.Ordinary; %finds R2 of fit
   coeff=linfit.Coefficients.Estimate; %reports coeffecients of fit
   lineslope(k)=coeff(2); %slope of the fit
   conf=coefCI(linfit); %confidence interval of the fit coeffecients
   lowerconf(k)=conf(2,1); %lower bound of confidence interval
   upperconf(k)=conf(2,2); %upper bound of confidence interval
   uppererror(k)=abs(lineslope(k)-lowerconf(k)); % lower bound error
   lowererrror(k)=abs(lineslope(k)-upperconf(k)); % lower bound error
     
end

for binfromend=0:k-1 %counts from the last bin to the first bin, 0 is the end/last bin, k-1 is the first bin, which is k-1 away from the last bin

 xfit =xg((1+stepsize*(k-binfromend)-stepsize):(stepsize*(k-binfromend)-stepsize+binwidth)); %fit region for time values for bin of interest
 yfit=yg((1+stepsize*(k-binfromend)-stepsize):(stepsize*(k-binfromend)-stepsize+binwidth));%fit region for radius values for bin of interes

   
 f = fit(xfit,yfit,'exp1');%fits to an expential model of form y=a*exp(bx) where a and b are coefficients
 resultmatrix=coeffvalues(f);%finds coeffeienct values 
 a1=resultmatrix(1); %a term in exponential model
 b1=resultmatrix(2);%b term in exponential model
 lambda=(b1*-3)^-1;%extensional relaxation time
 linf=fitlm(xfit,log(yfit)); % linear fit of radius and time data-EC region should be linear
 R_squaredbest=linf.Rsquared.Ordinary; %R squared of the fit- measure of goodness of fit

 
 
 if binfromend<k-1 %if the bin is not the first bin, then proceed
 for q=1:(k-binfromend-1) % looks at all the bins in front of the bin of interest.  (Ex. with 12 total bins, if you are 4 bins from the end (the 8th bin), then this looks at the 7 bins before this one)
     totalpointsn=(binwidth-stepsize)+((k-binfromend)-(k-binfromend-q)+1)*stepsize; % the total number of points to be binned
   xfitn =xg((1+stepsize*(k-binfromend-q)-stepsize):(1+stepsize*(k-binfromend-q)-stepsize+totalpointsn-1)); %identifies the points in the bin of interest and the q bins before the bin of interest for the time (x) values
 yfitn=yg((1+stepsize*(k-binfromend-q)-stepsize):(1+stepsize*(k-binfromend-q)-stepsize+totalpointsn-1));%identifies the points in the bin of interest and the q bins before the bin of interest for the radius (y) values
 
fn = fit(xfitn,yfitn,'exp1'); %fits to an expential model of form y=a*exp(bx) where a and b are coefficients
  resultmatrixn=coeffvalues(fn); %finds coeffeienct values 
 a1=resultmatrixn(1); %a term in exponential model
 b1=resultmatrixn(2);% b term in exponetial model
 linfn=fitlm(xfitn,log(yfitn));% linear fit of ln(radius) and time data-EC region should have a constant slope
 R_squaredbestn=linfn.Rsquared.Ordinary; % saves the R^2 of the lin fit
 ymodeln=a1.*exp(b1.*xfitn); %puts the coeffienct values in the model
 yresidualn=yfitn(end-stepsize:end)-ymodeln(end-stepsize:end);%compares the model with the data points
 yres_av(q)=mean(yresidualn);%averages the residuals
end
 locatebin=abs(yres_av)<tol; %checks if the average of the residuals is below the tolerance
 locatebin=find(locatebin==0);%locates the bins where resiudal is above tolerance
 if isempty(locatebin)==1 %if there are no bins above the tolerance
     binstart=1; % the bin start should be 1
 else
     binstart=k-binfromend-(locatebin(1)-1);% You have k bins total, locatebin-1 shows you the last bin to include that should still be within tol
 end
 end
binend=k-binfromend; % sets the bin end as the bin of interest defined at the beginning of the curent for loop
if binfromend>0 % if the bin from the end is not the last bin, then check all the bins following to see if they are within tolerance to include in the fit
    for q=1:binfromend %from 1 bin away from the last bin to the bin from the end
        totalpointsn1=(binwidth-stepsize)+((k-binfromend+q)-binstart+1)*stepsize;% the total number of points to be binned
   xfitn1 =xg((1+stepsize*(binstart)-stepsize):(1+stepsize*(binstart)-stepsize+totalpointsn1-1));%identifies the points in the bin of interest and the q bins after the bin of interest for the time (x) values
 yfitn1=yg((1+stepsize*(binstart)-stepsize):(1+stepsize*(binstart)-stepsize+totalpointsn1-1));%identifies the points in the bin of interest and the q bins after the bin of interest for the radius (y) values
 
fn1 = fit(xfitn1,yfitn1,'exp1');%fits to an expential model of form y=a*exp(bx) where a and b are coefficients
resultmatrixn1=coeffvalues(fn1); %finds coeffeienct values 
a1=resultmatrixn1(1); %a term in exponential model
b1=resultmatrixn1(2);%b term in exponential model
linfn1=fitlm(xfitn1,log(yfitn1));% linear fit of ln(radius) and time data-EC region should have a constant slope
R_squaredbestn=linfn1.Rsquared.Ordinary;% saves the R^2 of the lin fit
ymodeln1=a1.*exp(b1.*xfitn1); %puts the coeffienct values in the model
yresidualn1=yfitn1(end-stepsize:end)-ymodeln1(end-stepsize:end);%compares the model with the data points
yres_av1(q)=mean(yresidualn1); %averages the residuals
    end
    locatebin1=abs(yres_av1)<tol; %checks if the average of the residuals is below the tolerance
 locatebin1=find(locatebin1==0);%locates the bins where resiudal is above tolerance
 if isempty(locatebin1)==1%if there are no bins above the tolerance
     binend=k;%then the bin end should be the last bin, bin k
 else 
     binend=k-binfromend+locatebin1(1)-1;%else, the last bin should be the last bin where the fit matches within the tolerance
 end
 
end
binrange(binfromend+1,1)=binstart; %saves the bin start for the for loop centered around the bin of interest, binfromend
binrange(binfromend+1,2)=binend;%saves the bin end for the for loop centered around the bin of interest, binfromend
bindiff(binfromend+1)=abs(binrange(binfromend+1,1)-binrange(binfromend+1,2)); %finds how many bins are involved in the fit for the bin of interest
end
maxbindiff=find(bindiff==max(bindiff)); % finds the fit that uses the greatest number of bins
binstart=binrange(maxbindiff(1),1);%sets as the binstart (from fit that used greatest number of bins)
binend=binrange(maxbindiff(1),2);%sets as the binend (from fit that used greatest number of bins)
totalpoints=(binwidth-stepsize)+(binend-binstart+1)*stepsize; %fits over all the points in the full range
 xfit =xg((1+stepsize*binstart-stepsize):(1+stepsize*binstart-stepsize+totalpoints-1));%identifies x range (time range) of the fit
 yfit=yg(((1+stepsize*binstart-stepsize):(1+stepsize*binstart-stepsize+totalpoints-1)));%identifies y range (time range) of the fit
 
 f = fit(xfit,yfit,'exp1');%fits to an expential model of form y=a*exp(bx) where a and b are coefficients
 resultmatrix=coeffvalues(f);%finds coeffeienct values 
 a=resultmatrix(1); %a term in exponential model
 b=resultmatrix(2);%b term in exponential model
 lambda=(b*-3)^-1;%extensional relaxation time
 linf=fitlm(xfit,log(yfit));% linear fit of ln(radius) and time data-EC region should have a constant slope
 R_squaredbest=linf.Rsquared.Ordinary;% saves the R^2 of the lin fit
 coeff=linf.Coefficients.Estimate;%finds coeff
 linfslop=coeff(2);%model slope
 ymodel=a.*exp(b.*xfit);%model fit
 yresidual=yfit-ymodel;%finds residual

if binend==binstart %warns user if the fit was only on 1 bin
    msgbox('Caution:Fitting done over 1 bin only. May need to pick a different start point for binfromend or decrease your binwidth.');
end

% %Figures

figure (1) %a graph of the bins with slope- helpful to determine if binning is sensible- if user can change bin number by a few and fitting region dramatically changes, too few bins might be used
hold on
errorbar(linspace(1,binstart-1,binstart-1),lineslope(1:binstart-1),uppererror(1:binstart-1),lowererrror(1:binstart-1),'o','color','b','MarkerSize',10)
errorbar(linspace(binstart,binend,(binend-binstart+1)),lineslope(binstart:binend),uppererror(binstart:binend),lowererrror(binstart:binend),'o','color','r','MarkerSize',10)
if binend~=k
errorbar(linspace(binend+1,k,(k-binend-1+1)),lineslope(binend+1:k),uppererror(binend+1:k),lowererrror(binend+1:k),'o','color','b','MarkerSize',10)
end
xlabel('Time step','Interpreter','latex','FontSize',30)
ylabel('Slope','Interpreter','latex','FontSize',30)
title('Slope with confidence intervals','Interpreter','latex','FontSize',30)
set(gca,'FontSize',20)
hold off


figure (2) %shows the exponential decay region (EC regime) on a semilog scale. This should look linear, but it it doesn't, try changing binsize, tolerance or consider if you really have an EC region
plot(xfit,yfit,'MarkerSize',10)
set(gca,'YScale','log','FontSize',20);
% curtick = get(gca, 'XTick');
% set(gca, 'XTickLabel', cellstr(num2str(curtick(:))))
xlabel('time [s]','Interpreter','latex','FontSize',30)
ylabel('R [mm]','Interpreter','latex','FontSize',30);
title('$Exponential Decay, log(y)$','Interpreter','latex','FontSize',30)



figure(3) %shows the fit vs the data. Vertical lines are the bins.  Useful to see if need to lower tolerance- maybe be needed if the end or begginging slope bins does not match well with the rest of the fit

 axes('XScale', 'linear', 'YScale', 'log')

% legend('10mM CTAB/10mM NaSal')
title('Exponential Decay with Fitting function','Interpreter','latex','FontSize',30)
hold on
line([xfit(1),xfit(1)],[0.0001,1],'color','k','LineWidth',0.5)
for kk=1:2:length(xfit)/stepsize
   line([xfit(stepsize*kk),xfit(stepsize*kk)],[0.0001,1],'color','k','LineWidth',0.5) 
end
for kk=2:2:length(xfit)/stepsize
   line([xfit(stepsize*kk),xfit(stepsize*kk)],[0.0001,1],'color','k','LineWidth',0.5) 
end
  plot(f,xfit,yfit)
hold off
ylim([min(yfit) max(yfit)])
xlabel('time [s]','Interpreter','latex','FontSize',30)
ylabel('R/R$_{o}$','Interpreter','latex','FontSize',30');
set(gca,'FontSize',20)
yline=feval(f,xfit)


figure (4) %plots residuals- if the residuals are strongly trending away from the dotted line, may need to adjust tol
plotResiduals(linf,'probability','MarkerSize',10)
xlabel('Residuals','Interpreter','latex','FontSize',20);
ylabel('Probability','Interpreter','latex','FontSize',20);
title('Residuals,Probability Fit','Interpreter','latex','FontSize',20)
%legend('y(x)', sprintf('R^2 = %0.3f',R))
set(gca,'FontSize',30)


 
figure (5) %plots the radius by time over the total range (not just EC)
 plot(xtot,ytot)
set(gca,'YScale','log','Fontsize',20);
xlabel('time [s]','Interpreter','latex','FontSize',30)
ylabel('R/R$_{o}$','Interpreter','latex','FontSize',30)
title('Total data, y-axis log','Interpreter','latex','FontSize',30)


figure (6) %vizualization of residuals- if residuals trend strongly in one direction above or below red line, may need to lower tol
plot(xfit,yresidual)
title('Residuals','Interpreter','latex','FontSize',30)
xlabel('time [s]','Interpreter','latex','FontSize',30)
ylabel('Residual','Interpreter','latex','FontSize',30)
hold on

 line([0,xfit(end)], [0,0],'color','r')
 
line([xfit(1),xfit(1)],[0.0001,1],'color','k')
for kk=1:2:length(xfit)/stepsize
   line([xfit(stepsize*kk),xfit(stepsize*kk)],[-1,1],'color','k') 
end
for kk=2:2:length(xfit)/stepsize
   line([xfit(stepsize*kk),xfit(stepsize*kk)],[-1,1],'color','k') 
end
hold off
 ylim([min(yresidual) max(yresidual)])
xlim([min(xfit) max(xfit)])
 set(gca,'FontSize',20)
 hold off

% ax1 = gca; % current axes
% ax1_pos = ax1.Position; % position of first axes
% ax2 = axes('Position',ax1_pos,...
%     'XAxisLocation','top',...
%     'YAxisLocation','right',...
%     'Color','none','XTickLabel',bin1:1:binend);
% ax2.XColor = 'r';
% ax2.YColor = 'none';



T = table(endtime,R_squaredbest,linfslop,lambda,binstart,binend,binwidth,stepsize,a,b) %outputs a table of the endtime of thinning, the R^2 of the EC fit, the slope of the fit, the lamba or relaxation time in s, the start bin of the fit, the end bin of the fit, the bin width, step size, a prefeactor of the fit and b coeffecient of the fit. This table may be useful to save as a record of the data analysis