%INTRODUCTION
%This code uses raw exported data (with minimal file manipulation) from
%the Mastersizer 3000 software (Malvern, England) to calculate the fractal
%dimension of a single measurement timepoint, idendified below as 'c'. 
%
%See Farner Budarz et al. Langmuir, 2017 for a more detailed description
%and example of implementation. 
%
%


clear all
%figure(1)
%clf
%figure(2)
%clf
%figure(3)
%clf
%figure(4)
%clf


%BEFORE RUNNING, Prep data from mastersizer------------------
%Delete first row (column titles), and second column (sample name) of
%mastersizer file and save as tab delimeted "yyyymmdd.txt"

%load data from mastersizer----------------------------------
%xxxxxxxx.txt data from mastersizer where x=yyyymmdd

data=load('20150702.txt'); 

q=data(:,110:172); %rows vary w/ number of runs, always 63 columns (detectors)
i=data(:,173:235);

qt=transpose(q); %always 63 rows (detectors), columns vary w/ number of runs
it=transpose(i);

maxq = find(qt(:,1)==0); %identify zero values for scattering angle - these correspond to retrodiffusion angles
logq = log10(qt(1:maxq(1,1)-1,:)); %log transform of data (take only datapoints up to first zero value)
logi = log10(it(1:maxq(1,1)-1,:)); %same here to keep sizes the same

%Linear Regression of data

c=64; %select column (timepoint from Mastersizer data record)


figure(1) %plot full curve of log(q) vs log(i)
plot(logq(:,c),logi(:,c),'k--')
title({'Log(I) vs Log(q)' 'Linear Regression (Solid Blue), Measured Values (Dashed Black)'})
xlabel('log(q)')
ylabel('log(I)')

%Calculate slope of linear power region in log(q) vs log(i) plot
%
%Starting at the maximum q value and working backwards, this section
%compares the local slope of a 4 point segment of the scattering curve to
%the slope of the full range between the 4 point segment and maximum q. The
%goal is to determine the q value at which the local slope deviates from
%the slope of the full range, indicating the end of the the linear power
%law region. 
j=size(logq,1); %size of vector for for loop
for n=1:j-4; %for starting row position, begin at 
logqn = logq(n:n+4,c); %select 4 point segment of curve
login = logi(n:n+4,c);

lreg = polyfit(logqn,login,1); %linear regression of 4 point segment
%b(n,1) = lreg(2);
m(n,1) = lreg(1); %1st column of m is slope of 4 point segment

m(n,2) = (logi(end,c)-logi(n,c))/(logq(end,c)-logq(n,c)); % 2nd column of m is two point slope between first (n) (4 point segment) and the last q value
end

for z=1:size(m,1);
    m(z,3)=100*((m(z,1)-m(z,2))/m(z,1)); %3rd column of m is % diff between two slopes
end

figure(2) %show the calculated slopes
plot(logq(1:size(m,1),c),m(:,1),'k-')
hold on
plot(logq(1:size(m,1),c),m(:,2),'b')
title('Segment slope (black), first-last point slope (blue)')
xlabel('log(q)')
ylabel('Slope')

figure(3) %plot the percent difference between the calculated slopes
y=zeros(size(m,1),1);
plot(logq(1:size(m,1),c),y,'k-')
hold on
plot(logq(1:size(m,1),c),m(:,3),'bl')
axis([min(logq(1:size(m,1),c)) max(logq(1:size(m,1),c))+1 -25 25])
title('% Difference Between Slopes')
xlabel('log(q)')
ylabel('% Difference')

tolerance = 12; %identify the percent difference in slopes below which the two slopes are considered to be the same

%Here we start to work backwards from the last point and look for where
%the 4-point slope first deviates from the first-last slope by more than
%tolerance (%)

index1 = find(abs(m(:,3))>tolerance); %find all values where the absolute %difference is greater than the tolerance

found=find(index1>1);%find the index of all values in index1 
index=max(found); %find the last index of found - this is the last point of index 1. 

startlogq = logq(index1(index(1,1),1),c); %identify the logq value (x value) to start the linear portion of curve
dflogq=logq(index1(index(1,1),1)+1:size(logq,1),c);
dflogi=logi(index1(index(1,1),1)+1:size(logq,1),c);
lreg2 = polyfit(dflogq,dflogi,1); %linear regression of linear region of curve
df(1,c)=abs(lreg2(1))
lreg2Y = polyval(lreg2,dflogq);

%find r^2 for linear regression of power law region
yresid = dflogi - lreg2Y;
SSresid = sum(yresid.^2);
SStotal = (length(dflogi)-1) * var(dflogi);
rsq = 1 - SSresid/SStotal


figure(1)
hold on
plot(dflogq,dflogi,'k--')
plot(dflogq,lreg2Y,'b')
set(gcf,'Color','white')

figure(4)
hold on
plot(dflogq,dflogi,'ko')
plot(dflogq,lreg2Y)
title('Linear Regression (Blue) compared to Measured Values (circles)')
xlabel('log(q)')
ylabel('log(I)')

