% [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)');