## Copyright (C) 2013-2016 Claudio Girardi, in3otd@qsl.net
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 2
## of the License, or (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

## usage: calkit_model(FILENAME, TYPE [, SAVE_GRAPHS])
##
## compute the calibration standars coefficients to be used
## with 8753/8510 type VNAs from the measured S-parameters
##
## FILENAME is a 1-port S-parameters file in Touchstone format
## containing the measured "standard" data
## TYPE is the standard type to fit ("OPEN", "SHORT" or "LOAD")
## SAVE_GRAPHS (true or false) selects wheter to save the computed graphs to a file (default: false)
##
## !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
## ! NOTE: this script is still a WORK IN PROGRESS    !
## ! to fit your particular "standard" implementation !
## ! you might need to tweak the lines in the source  !
## ! code marked with "TWEAK"                         !
#  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
##
## References:
##
## [1]   Agilent, "Specifying Calibration Standards for the Agilent 8510 Network Analyzer," Application Note 8510-5B
## [2]   Agilent, "Specifying Calibration Standards and Kits for Agilent Vector Network Analyzers," Application Note 1287-11
##
## Example:
##
##  calkit_model('measured_open.s1p', 'OPEN')

## Version: 0.31

## Changelog:
##   0.31: corrected short model S-parameters calculations - thanks to hendorog for spotting this issue
##   0.3 : change fileparts() usage to work with Octave 4.x; minor cosmetical changes to plots
##   0.2 : corrected standard model losses, reorganized code
##   0.1 : initial version


function calkit_model(infile, type, save_graphs = false)

  if (nargin < 2)
    usage();
    error("Missing parameters.");
  endif

  global Z0;
  global c;
  global std_type;
  global s_meas;
  global xmin;
  global xmax;

  std_type = type;

  #debug_on_warning(false);

  # reference impedance
  Z0 = 50; # ohm
  # speed of light in vacuum
  c = 299792458; # m/s

  # directory containing measurements data
  [basedir, fbasename, fext] = fileparts(infile);

  graph_title = fbasename
  graph_title(graph_title == "_") = " " # remove underscores
  graph_title = [graph_title ' - ' toupper(std_type)];

  # load measured S-parameters
  [freq, data] = SXPParse(infile);
  s_meas.freq = freq';
  s_meas.s11 = squeeze(data(1,1,:));
  omega = 2 * pi * s_meas.freq;

  # fit between 10 MHz and 6 GHz
  xmin = lookup(s_meas.freq, 50e6); # (TWEAK)
  xmax = lookup(s_meas.freq, 12e9);  # (TWEAK)

  # Z11=Zo*(1+S11)/(1-S11)
  Z11 = Z0 * (1+s_meas.s11)./(1-s_meas.s11);
  X11 = imag(Z11);
  R11 = real(Z11);

  Y11 = 1 ./ Z11;
  G11 = real(Y11);
  B11 = imag(Y11);

  fig_idx = 1;
  figure(fig_idx);

  switch tolower(std_type)
    case "open"
      # plot Rp and Cp
      loglog(s_meas.freq, abs(1 ./ G11), 'r;Rp;');
      hold on;
      loglog(s_meas.freq, abs(1e15*B11) ./ omega, 'b;Cp;');
      hold off;
      title([graph_title ' - Rp and Cp - measured'], 'Interpreter','tex');
      xlabel('frequency, Hz');
      ylabel('Rp, ohm / Cp, fF');
      grid on;
    case "short"
      # plot Rs and Ls
      loglog(s_meas.freq, R11 , 'r;Rs;');
      hold on;
      loglog(s_meas.freq, abs(1e9*X11) ./ omega, 'b;Ls;');
      hold off;
      title([graph_title ' - Rs and Ls - measured'], 'Interpreter','tex');
      xlabel('frequency, Hz');
      ylabel('Rs, ohm / Ls, nH');
      grid on;
    case "load"
      # plot Rs and Xs
      loglog(s_meas.freq, R11 , 'r;Rs;');
      hold on;
      loglog(s_meas.freq, X11, 'b;Xs;');
      hold off;
      title([graph_title ' - Rs and Xs - measured'], 'Interpreter','tex');
      xlabel('frequency, Hz');
      ylabel('Rs, ohm / Xs, ohm');
      grid on;
    otherwise
      usage()
      error("Specified standard type (%s) not recognized.\n", std_type);
  endswitch

  ### optimization
  #
  # starting values for the offset line
  offs_delay = 35; # ps (as in the standard definition table) (TWEAK)
  offs_loss = 0.7; # Gohm/s (as entered in the standard definition table) (TWEAK)
  offs_Z0 = 50; # ohm (TWEAK)

  switch tolower(std_type)
    case "open"
      # starting values for the actual open capacitance
      C0 = 50e0; # fF (TWEAK)
      C1 = -310e0; # 1e-27 F/Hz (TWEAK)
      C2 = 23e0; # 1e-36 F/Hz^2 (TWEAK)
      C3 = -159e-3; # 1e-45 F/Hz^3 (TWEAK)

      # starting point parameters
      start_params = [offs_delay, offs_loss, offs_Z0, C0, C1, C2, C3];

      opt.lower_bounds = [10, 0.1, 30, -1e4, -1e4, -1e4, -1e4]; # (TWEAK)
      #opt.lower_bounds = 1e-3 .* abs(start_params);
      #opt.lower_bounds = 0.0 * ones(size(start_params));

      opt.upper_bounds = [50, 50, 60, 1e4, 1e4, 1e4, 1e4]; # (TWEAK)
      #opt.upper_bounds = 1e3 .* abs(start_params);
      #opt.upper_bounds = +inf * ones(size(start_params));

    case "short"
      # for 8753 VNAs the short has no inductance parameters
      #  (TWEAK if your VNA supports the inductance parameters for the short model)
      L0 = 0.0;
      L1 = 0.0;
      L2 = 0.0;
      L3 = 0.0;

      # starting point parameters
      start_params = [offs_delay, offs_loss, offs_Z0, L0, L1, L2, L3];

      #  (TWEAK if your VNA supports the inductance parameters for the short model)
      opt.lower_bounds = [20, 0.1, 10, 0.0, 0.0, 0.0, 0.0];
      opt.upper_bounds = [50, 5, 100, 0.0, 0.0, 0.0, 0.0];

    case "load"
      # starting point parameters
      # load has only the offset line parameters
      start_params = [offs_delay, offs_loss, offs_Z0];

      opt.lower_bounds = [20, 0.0, 10]; # (TWEAK)
      opt.upper_bounds = [100, 20, 100]; # (TWEAK)

    otherwise
      usage()
      error("Specified standard type (%s) not recognized.\n", std_type);
  endswitch

  # TWEAK (the default algorithms should work well, but you may want to experiment...)
  # global algos
  #opt.algorithm = NLOPT_GN_DIRECT_L;
  #opt.algorithm = NLOPT_GN_ORIG_DIRECT_L;
  #opt.algorithm = NLOPT_GN_DIRECT_L_RAND;
  opt.algorithm = NLOPT_GN_CRS2_LM;
  #
  #opt.algorithm = NLOPT_G_MLSL_LDS;
  #nlopt_opt_set_local_optimizer(opt, local_opt);
  #
  #opt.algorithm = NLOPT_GD_STOGO_RAND;
  #opt.algorithm = NLOPT_GN_ISRES;
  # local algos
  #opt.algorithm = NLOPT_LN_COBYLA;
  #opt.algorithm = NLOPT_LN_BOBYQA;
  # PRAXIS not well suited for bound constraints
  #opt.algorithm = NLOPT_LN_PRAXIS;
  #opt.algorithm = NLOPT_LN_SBPLX;

  opt.min_objective = @nlopt_sfunc;
  opt.xtol_rel = 1e-4;
  opt.verbose = 1;
  opt.maxeval = 100000;

  # TWEAK (2nd optimization is often not needed)
  [opt_params, fmin, retcode] = nlopt_optimize(opt, start_params);
  [ret_name, ret_desc] = NLopt_err(retcode);
  printf('1st optimization ended with status: %s (%s)\n', ret_desc, ret_name);

  # polish with a local algorithm
  opt.algorithm = NLOPT_LN_BOBYQA;
  #opt_params = start_params;
  [opt_params, fmin, retcode] = nlopt_optimize(opt, opt_params);
  [ret_name, ret_desc] = NLopt_err(retcode);
  printf('2nd optimization ended with status: %s (%s)\n', ret_desc, ret_name);

  # print optimization results
  printf("\n");
  printf("Optimization results:\n");
  printf("offs_delay = %f ps\n", opt_params(1));
  printf("offs_loss = %f Gohm/s\n", opt_params(2));
  printf("offs_Z0 = %f ohm\n", opt_params(3));

  switch tolower(std_type)
    case "open"
      printf("C0 = %f * 1e-15 F\n", opt_params(4));
      printf("C1 = %f * 1e-27 F/Hz\n", opt_params(5));
      printf("C2 = %f * 1e-36 F/Hz^2\n", opt_params(6));
      printf("C3 = %f * 1e-45 F/Hz^3\n", opt_params(7));
    case "short"
      printf("L0 = %f * 1e-12 H\n", opt_params(4));
      printf("L1 = %f * 1e-24 H/Hz\n", opt_params(5));
      printf("L2 = %f * 1e-33 H/Hz^2\n", opt_params(6));
      printf("L3 = %f * 1e-42 H/Hz^3\n", opt_params(7));
    case "load"
      1; # nothing
    otherwise
      usage()
      error("Specified standard type (%s) not recognized.\n", std_type);
  endswitch

  printf("\n");
  printf("Residual RMS error : %f dB\n", 20*log10(fmin))

  # compute model reflection coefficient
  s_mod.s11 = std_model_s_params(opt_params);

  sdiff = s_meas.s11 - s_mod.s11;

  fig_idx = fig_idx + 1;
  # plot S11
  figure(fig_idx);
  semilogx(s_meas.freq, 20*log10(abs(s_meas.s11)), 'r;measured;');
  hold on;
  semilogx(s_meas.freq, 20*log10(abs(s_mod.s11)), 'b;model;');
  hold off;
  title([graph_title ' - S_{11}'], 'Interpreter','tex');
  xlabel('frequency, Hz');
  ylabel('S_{11}, dB');
  axis([300e3 6e9 -1 1],"autoy")
  legend("boxoff");
  grid on;
  if (save_graphs)
    print(strcat(fbasename, '_S11_dB_meas_mod.png'), '-dpng', '-F:12');
  endif

  fig_idx = fig_idx + 1;
  # polar plot (sort of)
  figure(fig_idx);
  # limit plot to xmax
  plot(real(s_meas.s11(1:xmax)), imag(s_meas.s11(1:xmax)), 'r+;measured;');
  hold on;
  plot(real(s_mod.s11(1:xmax)), imag(s_mod.s11(1:xmax)), 'b+;model;');
  hold off;
  title([graph_title ' - S_{11} - measured vs. model'], 'Interpreter','tex');
  xlabel('real(S_{11})');
  ylabel('imag(S_{11})');
  # force square aspect ratio
  axis([-1 1 -1 1],"square")
  legend("boxoff");
  switch tolower(std_type)
    case "open"
      1; # nothing
    case "short"
      legend("location", "southwest");
    case "load"
      1; # nothing
  endswitch
  grid on;
  if (save_graphs)
    print(strcat(fbasename, '_S11_meas_mod.png'), '-dpng', '-F:12');
  endif

  fig_idx = fig_idx + 1;
  # S-parameters delta
  figure(fig_idx);
  #plot(real(sdiff), imag(sdiff), 'r;S11 residual;');
  semilogx(s_meas.freq, 20*log10(abs(sdiff)), 'r;residual;');
  title([graph_title ' - S_{11} modeling residual'],  'Interpreter','tex');
  xlabel('frequency, Hz');
  ylabel('\Delta S_{11}, dB');
  axis([300e3 6e9 -1 1],"autoy")
  legend("boxoff");
  grid on;
  if (save_graphs)
    print(strcat(fbasename, '_delta_S11_meas_mod.png'), '-dpng', '-F:12');
  endif
endfunction

function s11 = std_model_s_params(params)
  # compute the standard input reflection coefficient
  #   according to (1.4) in [2]
  global std_type;
  global Z0; # reference impedance
  global c; # speed of light
  global s_meas;

  freq = s_meas.freq;
  omega = 2 * pi * freq;

  offs_delay = params(1) * 1e-12; # convert from ps to s
  offs_loss = params(2) * 1e9; # convert from Gohm/s to ohm/s
  offs_Z0 = params(3);

  switch tolower(std_type)
    case "open"
      C0 = params(4);
      C1 = params(5);
      C2 = params(6);
      C3 = params(7);
      # from eq (1.13) in [2], scaling as in (1.15)
      C = C0 * 1e-15 + ((((C3*1e-45.*freq) + C2*1e-36).*freq) + C1*1e-27).* freq;
      # from eq (1.1) and (1.13) in [2]
      gamma_t = ((1 - j * omega .* C * Z0) ./ (1 + j * omega .* C * Z0));
    case "short"
      L0 = params(4);
      L1 = params(5);
      L2 = params(6);
      L3 = params(7);
      # from eq (1.12) in [2], scaling as in (1.16)
      L = L0 * 1e-12 + ((((L3*1e-42.*freq) + L2*1e-33).*freq) + L1*1e-24).*freq;
      # from eq (1.1) and (1.12) in [2]
      gamma_t = (j * omega .* L - Z0) ./ (j * omega .* L + Z0);
    case "load"
      gamma_t = 0.0; # assume 50 ohm load
    otherwise
      print_usage();
      error("Specified standard type (%s) not recognized.\n", std_type);
  endswitch

  # from eq (1.10) in [2]
  alpha_l = sqrt(freq ./ 1e9) * (offs_loss * offs_delay) / (2 * offs_Z0);
  beta_l = omega * offs_delay + alpha_l;
  gamma_l = alpha_l + j*beta_l;
  Zc = offs_Z0 + (1-j) .* sqrt(freq ./ 1e9) .* (offs_loss ./ (2 * omega));

  # from eq (1.1) in [2]
  gamma_1 = (Zc - Z0) ./ (Zc + Z0); # offset line reflection coeff.

  exp_m2gamma_l = exp(-2 * gamma_l);

  # eq (1.4) in [2]
  s11 = (gamma_1 .* (1 - exp_m2gamma_l - gamma_1 .* gamma_t) + exp_m2gamma_l .* gamma_t) ./ (1 - gamma_1 .* (exp_m2gamma_l .* gamma_1 + gamma_t .* (1 - exp_m2gamma_l)));
endfunction;

# objective function for NLopt
# see http://ab-initio.mit.edu/wiki/index.php/NLopt_Matlab_Reference
function [ssse, gradient] = nlopt_sfunc(params)
  # fit to S-parameters
  # TWEAK (you can define your own measure of fitting error)
  global s_meas;
  global xmin;
  global xmax;

  if (nargout > 1)
    gradient = 0;
  end

  # get the standard model S-parameters
  s_mod.s11 = std_model_s_params(params);

  # compute difference between model and measured data
  # over the predefined frequency range
  s_mod_1 = s_mod.s11(xmin:xmax);
  s_meas_1 = s_meas.s11(xmin:xmax);
  sdiff_1 = s_mod_1 - s_meas_1;
  # root-mean-square error (RMSE)
  # (other error measures could also be used)
  ssse = sqrt((sdiff_1' * sdiff_1) / length(sdiff_1));
endfunction

# NLopt return values description
# from http://ab-initio.mit.edu/wiki/index.php/NLopt_Reference#Return_values
function [ret_name, ret_desc] = NLopt_err(ret_code)
  # Successful termination (positive return values)
  switch (ret_code)
    case 1
      ret_name = 'NLOPT_SUCCESS';
      ret_desc = 'Generic success return value.';
    case 2
      ret_name = 'NLOPT_STOPVAL_REACHED';
      ret_desc = 'Optimization stopped because stopval was reached.';
    case 3
      ret_name = 'NLOPT_FTOL_REACHED';
      ret_desc = 'Optimization stopped because ftol_rel or ftol_abs was reached.';
    case 4
      ret_name = 'NLOPT_XTOL_REACHED';
      ret_desc = 'Optimization stopped because xtol_rel or xtol_abs was reached.';
    case 5
      ret_name = 'NLOPT_MAXEVAL_REACHED';
      ret_desc = 'Optimization stopped because maxeval was reached.';
    case 6
      ret_name = 'NLOPT_MAXTIME_REACHED';
      ret_desc = 'Optimization stopped because maxtime was reached.';
    # Error codes (negative return values)
    case -1
      ret_name = 'NLOPT_FAILURE';
      ret_desc = 'Generic failure code.';
    case -2
      ret_name = 'NLOPT_INVALID_ARGS';
      ret_desc = 'Invalid arguments (e.g. lower bounds are bigger than upper bounds, an unknown algorithm was specified, etcetera).';
    case -3
      ret_name = 'NLOPT_OUT_OF_MEMORY';
      ret_desc = 'Ran out of memory.';
    case -4
      ret_name = 'NLOPT_ROUNDOFF_LIMITED';
      ret_desc = 'Halted because roundoff errors limited progress. (In this case, the optimization still typically returns a useful result.)';
    case -5
      ret_name = 'NLOPT_FORCED_STOP';
      ret_desc = 'Halted because of a forced termination: the user called nlopt_force_stop(opt) on the optimization''s nlopt_opt object opt from the user''s objective function or constraints.';
    otherwise # should not happen
      ret_name = 'UNKNOWN_ERROR';
      ret_desc = 'Unknown error.';
  endswitch
endfunction

function usage
  printf('usage: calkit_model(FILENAME, TYPE [, SAVE_GRAPHS])\n');
  printf('\n');
  printf('compute the calibration standars coefficients to be used\n');
  printf('with 8753/8510 type VNAs from the measured S-parameters\n');
  printf('\n');
  printf('FILENAME is a 1-port S-parameters file in Touchstone format\n');
  printf('containing the measured "standard" data\n');
  printf('TYPE is the standard type to fit ("OPEN", "SHORT" or "LOAD")\n');
  printf('SAVE_GRAPHS (true or false) selects wheter to save the computed graphs to a file (default: false)\n');
  printf('\n');
  printf('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n');
  printf('! NOTE: this script is still a WORK IN PROGRESS    !\n');
  printf('! to fit your particular "standard" implementation !\n');
  printf('! you might need to tweak the lines in the source  !\n');
  printf('! code marked with "TWEAK"                         !\n');
  printf('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n');
  printf('\n');
endfunction
