%INTRODUCTION
%This code uses raw exported data (with minimal file manipulation) from
%the Mastersizer 3000 software (Malvern, England) to calculate the fractal
%dimension of each measurement timepoint in the exported file.
%
%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
%figure(5)
%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('20151008.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
%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.
for c=1:size(logq,2) %set up for loop for each column in logq (timepoint from Mastersizer data record)
figure(1) %plot logI vs logq
hold on % will plot each run
plot(logq(:,c),logi(:,c))
title('Log(I) vs Log(q)')
j=size(logq,1); %size of vector for for loop
for n=1:j-4; %for starting row position
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) and the last points
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)')
%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')
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);
size(index1);
index1(1,1);
end
Tdf=transpose(df)
%figure(4)
%hold on
%plot(dflogq,dflogi,'bo')
%plot(dflogq,lreg2Y)
figure(5) %Plot fractal dimension for each measurement number in original file
plot(df)
xlabel('measurement #')
ylabel('Df')
%manually save Df data into excel spreadsheet