% Maxine Jonas - BEAM at MIT - February 2005
% -------------------------------------------------------------------------
% Establishes a bead's mean squared displacement from its trajectory.
% -------------------------------------------------------------------------
% Modified by Heejin Choi for BE309 class, Oct 23, 2006
% 
% Format of File : txt file without any heading
%                   First column ==> x position
%                   Second column ==> y position
%                   Unit : um
% Size = # of data point
% SamplingFreq = 30 [Frame/sec] of CCD camera
% Marker = 'o'

%function GetMSD_BE309 (File, Size, SamplingFreq)

close all, clear all

File = 'subtracted_traj_sum30.txt';
Size = 178;
SamplingFreq = 1;
EntireMtx = dlmread (File,'\t'); 
'EntireMtx filled up.'


TauVector = [[1:1:19] double(int32(logspace(log10(20),log10(Size), 46)))];
% and make sure that length (TauVector) is a power of 2 for Matlab
L = length(TauVector); PowerStop = floor(log2(L)); clear L
TauVector = TauVector (1, 1:2^PowerStop); L = length (TauVector);
% 'TauVector created.'

% Initialization
SquaredX = zeros (Size, L); SquaredY = zeros (Size, L); 
SquaredDisplacement = zeros (Size, L);
MSD = zeros (1, L); 
MSDnm = zeros (1, L);  
MSDum = zeros (1, L);
MSDtoSave = zeros (1,L);
'Matrices initialized.'

for i = 1 : L
    Tau = TauVector (1,i);
    for j = 1 : Size - Tau
        SquaredX (j,i) = (EntireMtx(j,1) - EntireMtx(j + Tau,1)).^2;
        SquaredY (j,i) = (EntireMtx(j,2) - EntireMtx(j + Tau,2)).^2;
        SquaredDisplacement (j,i) = SquaredX (j,i) + SquaredY (j,i); % in um^2.
    end
    MSD(1, i) = 1e-12 * mean (real(SquaredDisplacement(1 : Size - Tau, i))); % in m^2.
end
MSDnm (1,:) = 1e18 .* MSD(1,:); % in nm^2       
'MSD computed.'

MSDtoSave (2,:) = MSDnm (1,:);
MSDtoSave (1,:) = TauVector (1,:) ./ SamplingFreq;
y = [MSDtoSave(1,:) ; MSDtoSave(2,:)];
fid = fopen('MSD.txt','w');
fprintf(fid,'%6.3f \t %6.3f \n',y);
fclose(fid)
'MSD saved in MSD.txt'

figure(1)
loglog (TauVector ./ SamplingFreq, MSDnm, 'o-');
title (['Mean squared displacement - ', File]); ylabel ('MSD (nm^2)'); xlabel ('Tau (s)');

figure(2)
plot(EntireMtx(:,1), EntireMtx(:,2),'o-')  
title(' Endocytosed Bead Position Tracking ')
xlabel(' X position [\mum] ')
ylabel(' Y position [\mum] ')

% MSDum (1,:) = 1e12 .* MSD(1,:);
% figure(1)
% loglog (TauVector ./ SamplingFreq, MSDum, Marker);
% title (['Mean squared displacement - ', File]); ylabel ('MSD (microns^2)'); xlabel ('Tau (s)');