% [oscfit, coeff] = AnalyzeCoordinates(x, calibration, parameters)
%
% Calculates oscillator fit and parameters from droplet coordinates, and 
% plots the coordinates with the fit. The plot can be used to quickly
% assess the quality of the fit.
%
% x = Droplet coordinates recorded at 1000 fps (mm)
% calibration = How many millimeters one pixel is.
% parameters = Structure containing starting parameters for fitting.
%
% oscfit = General harmonic oscillator model fitted to data.
% coeff = Structure containing the coefficients calculated from the fit:
%
% coeff.
%    k = Magnetic spring constant (N/m)
%    Fmu = CAH-force (N)
%    beta = Viscous damping coefficient (Ns/m)
%
% Success of the fitting depends on the starting value of each fitting
% parameter. Each of the parameters has a default value. If the fitting
% works with those, you don't need to define any of them. If there are
% problems, you can define any or all of them (e.g. "parameters.T=300". 
% Default values are showed in parentheses.
%
% parameters.
%    maxamplitude = Maximum oscillation amplitude used in the fitting, mm 
%         (10). Harmonic oscillator model only works when the magnetic
%         field is harmonic. This is only true near the magnet axis. As a
%         result, only oscillations with sufficiently small amplitude can
%         be used. You may increase this value if there are not enough
%         oscillations for good fitting (i.e. you get weird results)
%         However, you will also increase the error in the results at the 
%         same time.
%    C = Equilibrium point (i.e. position of the magnet axis), mm (0). 
%         You probably don't need to touch this.
%    T = Oscillation period, ms (200). First thing to change if fitting
%         fails. You can easily estimate the period by plotting the
%         droplet coordinates.
%    a_0 = Max peak oscillation amplitude, mm (10). Easy to estimate by
%          plotting droplet coordinates.
%    a = Coefficient relating to CAH-force (3e-2). Might need tuning.
%    g = Coefficient relating to viscous force(5e-1). Might need tuning.
%    x_0 = Offset (i.e. phase) (0) ms. You probably don't need to touch 
%          this.
%


function [oscfit, coeff] = AnalyzeCoordinates(varargin)
% First argument is the coordinates.
xpixels=varargin{1}; 

% Second argument is the calibration.
calibration=varargin{2};

% Convert pixels to millimeters.
x=xpixels*calibration;

% If some of the parameters are defined, use those values. Otherwise use
% the defaults. 
if length(varargin) > 2
    f=varargin{3};  % f stands for "fitting"
    if ~isfield(f,'maxamplitude')
        f.maxamplitude=10;
    end
    if ~isfield(f,'C')
        f.C=0;
    end
    if ~isfield(f,'T')
        f.T=200;
    end
    if ~isfield(f,'a_0')
        f.a_0=10;
    end
    if ~isfield(f,'a')
        f.a=5e-2;
    end
    if ~isfield(f,'g')
        f.g=3e-1;
    end
    if ~isfield(f,'x_0')
        f.x_0=0;
    end
else % If no parameters were defined, use defaults for all of them.
    f.maxamplitude=10;
    f.C=0;
    f.T=200;
    f.a_0=10;
    f.a=5e-2;
    f.g=3e-1;
    f.x_0=0;
end

cursor=1; % Helper variable for finding the min and max amplitude

% Searches the maximum amplitude (don't worry about details) 
while 1   
    high=max(x(cursor:end)); % Find the maximum droplet coordinate
    thigh=find(x(cursor:end)==high)+cursor-1; % Find the corresponding time
    low=min(x(cursor:end));  % Find the minimum droplet coordinate
    tlow=find(x(cursor:end)==low)+cursor-1; % Find the corresponding time
    
    % If the difference between max and min amplitude is lower than
    % maxamplitude, we have found our starting point.
    if high-low <= f.maxamplitude 
        cursor=min([thigh,tlow]);
        break;
    else
        cursor=max([thigh,tlow]); % Otherwise we move to next oscillation.
    end
end
first=cursor;

disp('Starting point resolved');

% When the droplet is almost stationary, dynamic wetting effects might
% invalidate the harmonic model. We discard small oscillations just in
% case.
minamplitude = 1;

% Same as previous loop, except we look for the minimum value.
while 1 
    high=max(x(cursor:end));
    xhigh=find(x(cursor:end)==high)+cursor-1;
    low=min(x(cursor:end));
    xlow=find(x(cursor:end)==low)+cursor-1;
    if high-low <= minamplitude
        cursor=min([xhigh,xlow]);
        break;
    else
        cursor=max([xhigh,xlow]);
    end
end
last=cursor;
x=x(first:last);

% Starting point for fitting needs to be larger than mean value. If this is
% not the case, the script simply inverts the values.
if (x(1) < mean(x)) 
    x=-x;  
end   
disp('Ending point resolved');

% Actual fitting. Solving x(t) from the equation of general harmonic
% oscillator leads to quite complicated function:
oscfittype=fittype('(-1)^(floor(2*(x-x_0)/T))*a_0*(sqrt(pi^2+g^2)/pi*((1+a)*exp(-g*(floor(2*(x-x_0)/T)))-2*a*(1-exp(-g*((floor(2*(x-x_0)/T))+1)))/(1-exp(-g)))*exp(-2*g/T*((x-x_0)-(floor(2*(x-x_0)/T))*T/2))*(cos(2*pi/T*((x-x_0)-(floor(2*(x-x_0)/T))*T/2))*pi/sqrt(pi^2+g^2) + sin(2*pi/T*((x-x_0)-(floor(2*(x-x_0)/T))*T/2))*g/sqrt(pi^2+g^2))+a) + C');
options=fitoptions('Method', 'NonlinearLeastSquares', 'StartPoint',[f.C,f.T,f.a,f.a_0,f.g,f.x_0],'MaxIter', 1000,'MaxFunEvals', 1000,'TolFun', 1e-16, 'TolX', 1e-16);
oscfit = fit(transpose(1:length(x)), transpose(x), oscfittype, options);
disp('Fitting done');

% Calculate parameters
coeff=CalcCoeff(oscfit);
disp('Parameters calculated');
disp(coeff);

% Plot coordinates and the fit.
figure;
plot(x,'k.');
hold on;
plot(oscfit,'r-');
xlabel('t (ms)');
ylabel('x (mm)');