% Maxine Jonas - BEAM at MIT - February 2005
% -------------------------------------------------------------------------
% From GetMSD_Record.m
% Following equations described in paper # 3-8, derives G*, G' & G'' from a
% bead's mean squared displacement (MSD), recorded by GetMSD_Record.m
% Mason TG (2000), Rheol. Acta 39: 371-378.
% -------------------------------------------------------------------------
% Modified by Heejin Choi for BE309 class, Oct 06, 2006
% 
% Format of File : txt file without any heading
%                   First column ==> x position
%                   Second column ==> y position
%                   Unit : um
% Size = # of data point
% SamplingFreq = 1 [Frame/sec] of CCD camera
% Marker = 'o'

% function GetMSD_Record (File, Size, SamplingFreq, Marker)

clear all;  close all;

% function GetGstar_Mason (File, Size, SamplingFreq, Marker, UseAvg)

File = 'MSD.txt';
Size = 10094;
SamplingFreq = 30;
Marker = 'o-';
MSDtoSave = dlmread (File,'\t');  % in nm^2
UseAvg = 0;

if UseAvg == 0 % Gstar of a single bead
    TauVector (1,:) = MSDtoSave (:,1)'; 
    MSD (1,:) = 1e-18 .* MSDtoSave (:,2)'; % in m^2
else % Gstar from the average MSD of several beads
    TauVector (1,:) = Avg (1, :);
    MSD (1,:) = 1e-18 .* Avg (2, :);
end
L = size (TauVector, 2)


OmegaVector = zeros (1, L);
LnMSD = zeros (1, L);
LnTau = zeros (1, L);
GlobalFit = zeros (1, L); 
Alpha = zeros (1, L);
GammaApp = zeros (1, L);
Denominator = zeros (1, L);
Gstar = zeros (1,L); Gstar_Record = zeros (2,L);
GpMason = zeros (1, L); 
GppMason = zeros (1, L); 
% 'Matrices initialized.'

%%% Now from MSD(tau) to G*(omega)
OmegaVector = 1 ./ TauVector; % from large (left) to small (right) Omega
kB = 1.3806503e-23; % Boltzmann constant in m2 kg s-2 K-1
T = 312; % Temperature in absolute degrees
Radius = 0.5e-6; % Bead radius in meters
Factor = kB * T / (pi * Radius);

LnMSD (1,:) = log (MSD(1,:)); % in m^2.
LnTau (1,:) = log (TauVector(1,:));
% Below, I've replaced GammaPiece by GammaApp, the approximation of the
% gamma function presented in Mason, Rheol Acta (2000), # 3-8.
n = 5;
PolyGlobalFit = polyfit (LnTau, LnMSD, n);
GlobalFit = polyval (PolyGlobalFit, LnTau);
% 'MSD fitted.'

if(1) 
   figure(3)
   loglog(TauVector, MSD*1e18,'o-');
   xlabel('\tau'); ylabel('MSD [nm^2]'); title('MSD');
end

for k = 1 : n
    Alpha (1,:) = Alpha(1,:) + k * PolyGlobalFit (1, n + 1 - k) .* LnTau (1,:) .^ (k-1);
end
% 'Alpha computed.' % from small Tau = large Omega (left) to large Tau (right).

figure(4)
semilogx( OmegaVector(1,:), Alpha(1,:),'o-'  )
xlabel('\omega'); ylabel('\alpha'); title('\alpha vs \omega');

GammaApp (1,:) = 0.457 .* (1 + Alpha(1,:)).^2 - 1.36 .* (1 + Alpha(1,:)) + 1.90;
Denominator (1,:) = MSD (1,:) .* GammaApp(1,:) .* i.^(-Alpha(1,:));

Gstar (1,:) = Factor ./ Denominator (1,:);
GpMason (1,:) = abs (Gstar(1,:)) .* cos ((pi/2) .* Alpha(1,:));
GppMason (1,:) = abs (Gstar(1,:)) .* sin ((pi/2) .* Alpha(1,:));
% Gp = real(Gstar) and Gpp = imag (Gstar) return the same results

% Record it all now :)
Gstar_Record (1,:) = OmegaVector (1,:); Gstar_Record (2,:) = abs(Gstar (1,:));
Gp_Record (1,:) = OmegaVector (1,:); Gp_Record (2,:) = GpMason (1,:);
Gpp_Record (1,:) = OmegaVector (1,:); Gpp_Record (2,:) = GppMason (1,:);


'G*, Gp, Gpp saved'

figure(1)
loglog(OmegaVector, GpMason, Marker)
title(' Storage Modulus G prime ');
xlabel('\omega'); ylabel('pascals');

figure(2)
loglog(OmegaVector, GppMason, Marker)
title(' Loss Modulus G double prime ');
xlabel('\omega'); ylabel('pascals');