% File 'cal-model_fit'
% 
% For GNU Octave
% --------------------------------------------------------------------------------------------------
% 
% Octave script that fits VNA calibration standards parameters for open, short, match and thru to 
% measured S-Parameter data, uses the Keysight calibration standard parametrization
% 
% The one-port calibration standards are modeled by a lossy transmission line in front of a 
% frequency dependent impedance (capacitance for open, inductance for short, resistance for match). 
% The frequency dependence of the reactance is modeled by a polynomial of order three for open and 
% short. The transmission line is characterized by its offset wave impedance, offset delay and 
% offset loss. The approximations being made are described in the Keysight AN quoted below.
% 
% Literature:
% Keysight AN "Specifying Calibration Standards and Kits for Keysight Vector Network Analyzers"
% 
% Author: Dr. Mario Hellmich <mario.hellmich@web.de>
% Date:   04.02.2019
% 
% Uses the Octave Forge 'optim' package 
%
% --------------------------------------------------------------------------------------------------
% 
% 


1;                                                                           				   		% Prevent Octave from thinking this is a function file
clear;                                                                          					% Delete all
clc;																								% Clear screen
clf;																								% Clear plot
pkg load 'optim'																					% Load optim package


% ----- Constants -----
Z_0 = 50;																							% System impedance


% ----- Measured Data -----
trace_open  = importdata( 'eric_open.s1p', ' ', 5 );                       							% Returns struct, where .data contains the columns of the Touchstone file
trace_short = importdata( 'eric_short.s1p', ' ', 5 );
trace_match = importdata( 'eric_match.s1p', ' ', 5 ); 
trace_thru  = importdata( 'eric_thru.s2p', ' ', 5 ); 

s11_meas_open  = trace_open.data(:,2)  + i*trace_open.data(:,3);
s11_meas_short = trace_short.data(:,2) + i*trace_short.data(:,3);
s11_meas_match = trace_match.data(:,2) + i*trace_match.data(:,3);
s21_meas_thru  = trace_thru.data(:,4)  + i*trace_thru.data(:,5);

freq_open      = trace_open.data(:,1);
freq_short     = trace_short.data(:,1);
freq_match     = trace_match.data(:,1);
freq_thru      = trace_thru.data(:,1);


% ----- Cal Kit Model Definition -----
% Parameters open:  p(1): C_0,      p(2): C_1,     p(3): C_2,      p(4): C_3,
% Parameters short: p(1): L_0,      p(2): L_1,     p(3): L_2,      p(4): L_3,
% Parameters match: p(1): R
% Line parameters:  q(1): offset delay,   q(2): offset loss,       offset impedance is set to system impedance

% Impedances at the far end of the transmission line from the VNA port
scale_open  = [1e-15, 1e-27, 1e-36, 1e-45];															% Parameter scale open
scale_short = [1e-12, 1e-24, 1e-33, 1e-42];															% Parameter scale short
react_open  = @(p, f) p(1)*scale_open(1)  + p(2)*scale_open(2)*f  + ...
                      p(3)*scale_open(3)*(f.^2)  + p(4)*scale_open(4)*(f.^3);						% Frequency dependent reactance in open
react_short = @(p, f) p(1)*scale_short(1) + p(2)*scale_short(2)*f + ...
                      p(3)*scale_short(3)*(f.^2) + p(4)*scale_short(4)*(f.^3);						% Frequency dependent reactance in short
Z_open  = @(p, f) ( 1 / (2*pi*i) ) * ( 1 ./ ( f .* react_open(p, f) ) );							% Reactance in open standard
Z_short = @(p, f) ( 2*pi*i ) * ( f .* react_short(p, f) );											% Reactance in short standard
Z_match = @(p, f) p(1);																				% Resistance in match standard

% Parameter scale line
scale_l = [1e-12, 1e9];

% Wave impedance of transmission line
Z_c = @(q, f) ( Z_0 + q(2)*scale_l(2) * ( 1 ./ (4*pi*sqrt(f*1e9)) ) ) - i * ( q(2)*scale_l(2) * ( 1 ./ (4*pi*sqrt(f*1e9)) ) );

% Transmission constant gamma = alpha + i*beta times line length l, i.e. gamma * l
gammal = @(q, f) ( ( (q(2)*scale_l(2)*q(1)*scale_l(1))/(2*Z_0) ) * sqrt(f/1e9) ) + ...
                 i * ( 2*pi*f*q(1)*scale_l(1) + ( ( (q(2)*scale_l(2)*q(1)*scale_l(1))/(2*Z_0) ) * sqrt(f/1e9) ) );

% Reflection factor of an impedance Z with respect to system impedance Z_0
refct_0 = @(Z) ( Z - Z_0 ) ./ ( Z + Z_0 );

% Reflection factor at VNA port when the lossy line is terminated by an impedance Z as a function of the line parameters
refct_1 = @(q, f, Z) ( refct_0(Z_c(q,f)) .* ( 1 - exp(-2*gammal(q,f)) - refct_0(Z_c(q,f)) .* refct_0(Z) ) + refct_0(Z) .* exp(-2*gammal(q,f)) ) ./ ...
                     ( 1- ( refct_0(Z_c(q,f)) .* ( refct_0(Z_c(q,f)) .* exp(-2*gammal(q,f)) + refct_0(Z) .* (1-exp(-2*gammal(q,f))) ) ) );

% S-parameters of cal kit model as a function of the parameter vector pvect  
s11_model_open  = @(pvect, f) refct_1(pvect(1:2), f, Z_open(pvect(3:6), f));
s11_model_short = @(pvect, f) refct_1(pvect(1:2), f, Z_short(pvect(3:6), f));
s11_model_match = @(pvect, f) refct_1(pvect(1:2), f, Z_match(pvect(3), f));
s21_model_thru  = @(pvect, f) exp( -gammal(pvect, f) );


% ----- Goodness of Fit Evaluation for Complex Data -----
function retval = corcoeff( y, model_y )
	sum_squares = sumsq( real(y) - mean(real(y)) ) + sumsq( imag(y) - mean(imag(y)) );
	residuals   = sumsq( model_y - y );
	retval      = sqrt( residuals / sum_squares );
endfunction


% ----- Initial Values for Regression -----
xin_open  = [8,      6,       -5,        -300,      23,     -0.1 ]';
xin_short = [8,      2,        2,        -100,       2,    -0.01 ]';
xin_match = [2,      3,       50]';
xin_thru  = [42,     3]';

options = optimset('TolX', 1e-5, 'TolFun', 1e-6, 'MaxIter', 1e5, 'MaxFunEvals', 1e8);				% Set precision for optimization

% ----- Parameter Regression for Open Standard -----
%residual_open = @(x) sumsq( s11_model_open(x, freq_open) - s11_meas_open );
residual_open = @(x) [ real( s11_model_open(x, freq_open) - s11_meas_open ); imag( s11_model_open(x, freq_open) - s11_meas_open ) ];
[x_open, ~, cvg_open, outp_open] = nonlin_residmin( residual_open, xin_open, options);
%[x_open, ~, cvg_open, outp_open] = fminsearch( residual_open, xin_open, options);


% ----- Parameter Regression for Short Standard -----
%residual_short = @(x) sumsq( s11_model_short(x, freq_short) - s11_meas_short );
residual_short = @(x) [ real( s11_model_short(x, freq_short) - s11_meas_short ); imag( s11_model_short(x, freq_short) - s11_meas_short ) ];
[x_short, ~, cvg_short, outp_short] = nonlin_residmin( residual_short, xin_short, options );
%[x_short, ~, cvg_short, outp_short] = fminsearch( residual_short, xin_short, options );


% ----- Parameter Regression for Match Standard -----
%residual_match = @(x) sumsq( s11_model_match(x, freq_match) - s11_meas_match );
residual_match = @(x) [ real( s11_model_match(x, freq_match) - s11_meas_match ); imag( s11_model_match(x, freq_match) - s11_meas_match ) ];
[x_match, ~, cvg_match, outp_match] = nonlin_residmin( residual_match, xin_match, options );
%[x_match, ~, cvg_match, outp_match] = fminsearch( residual_match, xin_match, options );


% ----- Parameter Regression for Trhu Standard -----
%residual_thru = @(x) sumsq( s21_model_thru(x, freq_thru) - s21_meas_thru );
residual_thru = @(x) [ real( s21_model_thru(x, freq_thru) - s21_meas_thru ); imag( s21_model_thru(x, freq_thru) - s21_meas_thru ) ];
[x_thru, ~, cvg_thru, outp_thru] = nonlin_residmin( residual_thru, xin_thru, options );
%[x_thru, ~, cvg_thru, outp_thru] = fminsearch( residual_thru, xin_thru, options );


% ----- Print Results -----
result_open    = x_open  .* [scale_l, scale_open]'
cvg_open
corcoeff_open  = corcoeff( s11_meas_open, s11_model_open(x_open, freq_open) )

result_short   = x_short .* [scale_l, scale_short]'
cvg_short
corcoeff_short = corcoeff( s11_meas_short, s11_model_short(x_short, freq_short) )

result_match   = x_match .* [scale_l, 1]'
cvg_match
corcoeff_match = corcoeff( s11_meas_match, s11_model_match(x_match, freq_match) )

result_thru    = x_thru  .* [scale_l]'
cvg_thru
corcoeff_thru  = corcoeff( s21_meas_thru, s21_model_thru(x_thru, freq_thru) )


% ----- Plot Results -----
plot( s11_meas_open );
hold on;
plot( s11_model_open(x_open, freq_open) );
hold on;
plot (s11_model_open(xin_open, freq_open) );
title('Open');
xlabel('Re');
ylabel('Im')
legend('measured', 'model', 'initial');
print -dpng open_regression-result.png

hold off;
clf;

plot( s11_meas_short );
hold on;
plot( s11_model_short(x_short, freq_short) );
hold on;
plot (s11_model_short(xin_short, freq_open) );
title('Short');
xlabel('Re');
ylabel('Im')
legend('measured', 'model', 'initial');
print -dpng short_regression-result.png

hold off;
clf;

plot( s11_meas_match );
hold on;
plot( s11_model_match(x_match, freq_match) );
hold on;
plot (s11_model_match(xin_match, freq_open) );
title('Match');
xlabel('Re');
ylabel('Im')
legend('measured', 'model', 'initial');
print -dpng match_regression-result.png

hold off;
clf;

plot( s21_meas_thru );
hold on;
plot( s21_model_thru(x_thru, freq_thru) );
hold on;
plot (s21_model_thru(xin_thru, freq_thru) );
title('Thru');
xlabel('Re');
ylabel('Im')
legend('measured', 'model', 'initial');
print -dpng thru_regression-result.png