clear all; clf;
% To enable use of 'fsolve':
pkg load optim

% This section creates a noisy curve to fit the function to: ========
% Here are the pre-noise exact values which fsolve will approximate:
R = 100;
C = 1e-9;
L = 1e-5;
% Store input f and Z as two separate vectors
% Column vector of f values
f = logspace(4, 8, 100)';

% Create a noisy complex curve in Vc
Zl = j*2*pi*f*L;
Zc = 1./(j*2*pi*f*C);
Vc = Zc ./ (R + Zl + Zc);
rn = real(Vc) .* (1 + .05*randn(size(f)));
in = imag(Vc) .* (1 + .05*randn(size(f)));
Vc = rn + j*in;
% ===================================================================

fit2 = @(x) fit(x,f,Vc);

% Initial solution guess (column vector)
vin = [100; 2e-9; 1.2e-5];
% Display, TolFun, and TolX options don't seem to have an effect:
options = optimset('Display', 'iter', 'TolFun', 1e-12);
[vout, fval, info, output]=fsolve(fit2, vin, options);

% Print out the results
% 'vout' is the solution vector
% 'fval' is the final function evaluation (can be used to evaluate how
%    good the fit is)
% 'info = 1' means solution converged, but 2 and 3 seem to be okay too
% 'output' lists some statistics about the process
vout
info
output

% Plot input data & best fit curve
Zl = j*2*pi*f*vout(3);
Zc = 1./(j*2*pi*f*vout(2));
fitVc = Zc ./ (vout(1) + Zl + Zc);
inVc = Vc;
subplot(2, 1, 1);
hold;
semilogx (f, abs(inVc), 'r.');
semilogx (f, abs(fitVc), 'b');
subplot(2, 1, 2);
hold;
semilogx (f, arg(inVc)*180/pi, 'r.');
semilogx (f, arg(fitVc)*180/pi, 'b');