%This script plots the complex viscosity of a fluid from the
%power spectrum of the thermal fluctuations of a particle trapped in the
%fluid. To demonstrate, a test data set of an optically trapped 10um 
%silica bead in water (trap stiffness 2.5 uN/m) is used and can be 
%accessed here: http://dx.doi.org/10.5525/gla.researchdata.827 
%This demonstration starts with performing an FFT on a time series. If 
%data is already in FFT or PSD form then the FFT step can be skipped. 
%Discussion of the theory used in this code can be found in Chapters 2 & 
%7 of the thesis which this file accompanies.
%a.mcglone.1@research.gla.ac.uk. 


oversampling=1;                 %oversampling factor: increase for more density in autocorrelation
dpoints=250;                    %No. of points in autocorrelation function
K = 2.5;                        %trap stiffness micronewtons/m
rad = 10;                       %bead radius micrometres 
scale = K/(6*pi*rad);           %complex viscosity scaling factor
B = dlmread('k25rad10.txt', '\t'); %time series of trapped bead in water
gt = [B(:,1) B(:,2)];
D = B(:,2)-mean(B(:,2));        %subtract mean
Fx = (abs(fft(D)));             %Fourier transform
PSRFx = Fx/length(Fx);
RFx = PSRFx(1:(length(PSRFx)/2)).^2; % create one sided PSD
torigx = gt(:,1)-gt(1,1) ;
dt = torigx(2);
maxfreq = 1/dt;
minfreq = 1/max(torigx);
freqfft = linspace(minfreq,maxfreq/2,length(RFx));
k = 1;
N = 2*length(RFx);
df= freqfft(2)-freqfft(1);
Fs = length(RFx)*df;
for k = 1:1:(N-1)/2
RFx(N+1-k)= conj(RFx(k+1)) ; % 1 <= k <= (N - 1)/2 create symmetrical PSD
end
tat = dt*(0:N-1);
A = N*ifft((RFx)); %Inverse fourier transform of symmetrical PSD (autocorrelation)
NA =(A(2:end)./trapz(RFx));  %Normalise autocorrelation with PSD area
a=NA(1:dpoints);             %Choose autocorrelation length 
torig = linspace(dt,dt*dpoints,dpoints).';
t = linspace(dt,dt*dpoints,(dpoints*oversampling)-1).';
g = spline(torig, a, t);
GData = zeros(dpoints*oversampling,3);
frange = linspace(1/(torig(dpoints)),maxfreq,dpoints*oversampling);
g0 = 1;    
eta = 0;  
for ww = 1:dpoints*oversampling      %this function performs fourier transform on the autocorrelation function and multiplies by i*omega
    w = frange(ww);
    fta = ((1i*w*g0 + (1-exp(-1i*w*t(1)))*(g(1)-g0)/t(1) + sum(diff(g)./diff(t).*(exp(-1i*w*(t(1:(size(t)-1))))-exp(-1i*w*(t(2:size(t))))))))./(-w.^2);
    GStar = fta .* 1i*w ;
    GData(ww,:) = [w real(GStar) imag(GStar)];
end
GTOT = complex(GData(:,2),(GData(:,3)));
GF = complex(GTOT./(1-GTOT)).*scale;                      
nF = (sqrt((real(GF).^2)+(imag(GF).^2))./frange.'); %Complex viscosity 

figure('units','normalized','outerposition',[0 0 1 1])
subplot(2,1,2)
loglog(frange,nF,'-b')
ylim([0.0001 0.01])
xlabel('\omega(Rad.s)');
ylabel('\eta*(Pa.s)');
subplot(2,1,1)
semilogx(torig,a,'o')
hold on
semilogx(t,g,'x')
xlabel('t /s');
ylabel('A(t)');