
%% Create a trajectory

D=0.1;          %(micron^2/s)
timestep=.05;     %TimeStep (s)
pixelsize=.1;   %(micron)
sz=128;
T=500;          %Time Frames
Dim=2;          %Number of Dimensions

%This creates a trajectory
displacements=sqrt(2*D*timestep)*randn(T,Dim);
x=cumsum(displacements(:,1));
y=cumsum(displacements(:,2));

figure  %This opens a new, blank figure
plot(x,y,'r','linewidth',2)
xlabel('X (micron)','fontsize',14)
ylabel('Y (micron)','fontsize',14)
axis equal

%mean(sum(displacements.^2,2))/(2*Dim*timestep)

%% The 'while' loop

%Calculate a Mean Square Displacement (MSD) from the trajectory
N=round(T/4);
delta=(1:N)';   %time interval in frames

Num=zeros(N,1);
MSD=zeros(N,1);

for tt=delta'    %make MSD upto 25% of total time
    ii=1;
    tmp_sum=0;
    tmp_N=0;
    while (ii+tt)<T
        d2=(x(ii+tt)-x(ii))^2+(y(ii+tt)-y(ii))^2;
        tmp_sum=tmp_sum+d2;
        tmp_N=tmp_N+1;
        ii=ii+tt;    %non-overlapping displacements
    end
    Num(tt)=tmp_N;
    MSD(tt)=tmp_sum/tmp_N;
end

deltaT=delta*timestep   %time inteval in seconds

figure  %This opens a new, blank figure
plot(deltaT,MSD,'r','linewidth',2)
xlabel('\Deltat (s)','fontsize',14)
ylabel('MSD (\mum^2)','fontsize',14)
axis equal

%% Fit the MSD

%The diffusion constant can be found be from the slope (a)
%D=a/(2*Dim)

%Choose the method fitting
% Use either 'polyfit','fit' or 'fminsearch'

fitmethod='polyfit'     %Only for polynomial
fitmethod='fit'         %Requires Curve Fitting Toolbox
fitmethod='fminsearch'  %Uses a minimization method

switch fitmethod
    
    case 'polyfit'  %only for polynomials
        
        n=1; %order of polynomial, n=2 means go up to x^2
        found_coef=polyfit(deltaT,MSD,n);
        fit_curve=polyval(found_coef,deltaT);
        D=found_coef(1)/(2*Dim);
        
    case 'fit'  % 'fit' can do arbitrary curve fit models with weighting
        
        %define weights
        W = double(Num./deltaT.^2)
        
        opts=fitoptions('Method','Nonlinear','StartPoint',[.01 10],'Weights',W,'Robust','on');
        %'Robust' ignores outliers in the fitting. Usually only used if you have
        %some reason for justifying ignoring outliers
        
        ftype=fittype('a*x+b','coefficients',{'a','b'},'options',opts);
        [fi, goodnessL,output]=fit(deltaT,MSD,ftype)
        a_found=fi.a;
        b_found=fi.b;
        fit_curve=a_found*deltaT+b_found;
        D=found_coef(1)/(2*Dim);
        
    case 'fminsearch'
        
        %define weights
        W = double(Num./deltaT.^2)
        
        minfunction=@FitFunc;   %The function to be minimized
        options = optimset('Display','iter');
        
        %setting an initial guess
        X0(1)=1;
        X0(2)=0;
        
        %fminsearch iterates to minimize 'fval'
        [params,fval]=fminsearch(minfunction,X0,options,deltaT,MSD); %everything after 'options'
        %are the other inputs to the function.
        [fval fit_curve]=FitFunc(params,deltaT,MSD);
        D=params(1)/(2*Dim);
        
    otherwise
        fprintf('No valid fit method given.\n')
end

fprintf('Found Diffusion Coefficient: %g micron^2/s.\n',D)

figure  %This opens a new, blank figure
plot(deltaT,MSD,'r','linewidth',2)
hold on
plot(deltaT,fit_curve,'g','linewidth',2)
xlabel('\Deltat (s)','fontsize',14)
ylabel('MSD (\mum^2)','fontsize',14)
legend('Data','Fit')
title(fitmethod,'fontsize',16)
axis equal

%% Stuctures, Classes and Objects

% A simple stucture can be a nice data container

SPT.x=x;
SPT.y=y;
SPT.timestep=timestep;
SPT

% A class can also define functions

%an 'object' is an instance of a 'class'

%see SPTclass

SPTobj=SPTclass(x,y,timestep)

SPTobj.x

SPTobj.PlotTrajectory








