MATLAB_Utilities/fit_Nyquist_plot.m at master · Campbell-Muscle-Lab/MATLAB_Utilities · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
function [best_p,fit_em,fit_vm] = fit_Nyquist_plot(varargin);

% Handle inputs
params = inputParser;
addRequired(params,'freq');
addRequired(params,'elastic_mod');
addRequired(params,'viscous_mod');
addOptional(params,'scaling_vector',[1e6 0.1 1e6 10 1e6 10]);
addOptional(params,'initial_p',[1 1 1 1 1 1]);
addOptional(params,'lower_bounds',zeros(1,6));
addOptional(params,'upper_bounds',inf*ones(1,6));
addOptional(params,'figure_number',0);
parse(params,varargin{:});
params = params.Results;

% Code
expt_data.freq = params.freq;
expt_data.elastic_mod = params.elastic_mod;
expt_data.viscous_mod = params.viscous_mod;

best_p = fminsearchbnd( ...
    @return_Nyquist_fit, ...
    params.initial_p, ...
    params.lower_bounds,params.upper_bounds, ...
    optimset('MaxFunEvals',15000), ...
    expt_data, ...
    params.scaling_vector, ...
    0)

Nyquist_values = return_Nyquist_values(best_p, expt_data.freq, params.scaling_vector);
fit_em = real(Nyquist_values);
fit_vm = imag(Nyquist_values);

plot_fit(params.figure_number, expt_data, [], best_p, params.scaling_vector);

% Store for output
best_p = best_p .* params.scaling_vector;

end

function Nyquist_values = return_Nyquist_values(p, freq, scaling_vector)

    x_data = 2*pi*1i*freq;

    p = p.*scaling_vector;

    Nyquist_values = ...
        p(1) * (x_data .^ p(2)) - ...
        p(3) * (x_data ./ (2*pi*p(4) + x_data)) + ...
        p(5) * (x_data ./ (2*pi*p(6) + x_data));
end

function error_value = return_Nyquist_fit(p, expt_data, scaling_vector, figure_display)

    % Calculate Nyquist_values
    Nyquist_values = return_Nyquist_values(p,expt_data.freq, scaling_vector);

    % Split into elastic and viscous components
    fit_data.elastic_mod = real(Nyquist_values);
    fit_data.viscous_mod = imag(Nyquist_values);

    % Calculate sum of squares
    error_value = ...
        sum((fit_data.elastic_mod - expt_data.elastic_mod).^2) + ...
        sum((fit_data.viscous_mod - expt_data.viscous_mod).^2);

    % Display if required
    if (figure_display)
        plot_fit(figure_display, expt_data, fit_data);
    end

end

function plot_fit(figure_display, expt_data, fit_data, p, scaling_vector)

    figure(figure_display);
    cla;
    hold on;
    plot(expt_data.elastic_mod, expt_data.viscous_mod,'k^');

    if (isempty(fit_data))
         % Calculate Nyquist_values
        Nyquist_values = return_Nyquist_values(p,expt_data.freq, scaling_vector);
        % Split into elastic and viscous components
        fit_data.elastic_mod = real(Nyquist_values);
        fit_data.viscous_mod = imag(Nyquist_values);
    end

    plot(fit_data.elastic_mod, fit_data.viscous_mod,'bo-');
%     for i=1:numel(fit_data.elastic_mod)
%         plot([fit_data.elastic_mod(i) expt_data.elastic_mod(i)], ...
%             [fit_data.viscous_mod(i) expt_data.viscous_mod(i)],'r-');
%     end
end