%phase_tracking.m
%[e, u, zr, n, pro, X, Yn, Yd, m1, m2] = phase_tracking(f, s11, s21, d)
% Written by Jacob Lewis May 2020.
% The purpose of this program is to take the S-parameters accross a
% spectrum of frequencies to calculate the C-parameters of a substrate,
% compensating for branches and accounting for every quadrant of
% metamaterial. Though all values returned are complex, an assumption of
% the substrate is that it is lossless and thus the values may not be
% accurate if they indicate a level of loss or gain (such as imaginary
% permittivity, permeability or negative impedance). All values returned
% will be arrays of type "complex double". All arrays given and returned
% should be the same size (except for d which is scalar). The S-parameters
% should be complex doubles and f and d should be real.
% 
% Given:
%   f   Frequency
%   s11 S11
%   s21 S21
%   d   Substrate thickness (scalar,real)       *OPTIONAL, Default: 0.01 *
% Return:
%   e   Permittivity
%   u   Permeability
%   n   Index of refraction
%   z   Impedance
%   pro Wave propagation < e^(j \gamma d) >
%   X   Used to compute egd <egd = X ± sqrt(X^2 - 1)>
%   m1  Sign branch
%   m2  Log branch

function [e, u, zr, n, pro, X, Yn, Yd, m1, m2] = phase_tracking(f, s11, s21, d)
    
    if nargin < 4
        d = 0.01;   % d = 1 [cm]
    end
    
    % Start with these intial conditions
    m1 = ones(size(f)); % Sign branch is positive
    m2 = m1*0; % Log branch is 0

    % Calculate C-parameters
    [~, ~, pro, X, ~] = CParam(f, s11, s21, d, m1, m2);

    % Find sign branches
    phi = angle(pro);   % Angle of wave propagation
    phin = abs(phi-pi/2);   % Normalize to pi/2 and take abs
    % Local peaks above 1.5 of phin are the sign branches
    [~, pks] = findpeaks(max(phin,1.5));    % Signal Processing Toolbox
    m1(pks) = -1;       % Assign the peaks to -1 as inflection points
    m1 = cumprod(m1);   % Carry changes throughout
    m1 = (m1 + 1).*(real(X) < 1) - 1;   % m1 = -1 for ENG/MNG

    % Recalculate Sign Branch Corrected C-parameters
    [~, ~, pro, ~, ~] = CParam(f, s11, s21, d, m1, m2);

    % Find log branches
    phi = angle(pro);   % Angle of wave propagation
    % Local peaks above 3 of phi are the log branches
    [~, pks] = findpeaks(max(phi,3));
    val = phi(pks - 1) > 0; % val is for determining the sign of the branch
    % ± branch count ['just after', 'at'] branch location
    m2(pks + val) = 2*val - 1;
    m2 = cumsum(m2);        % Carry changes throughout
    
    % Check if DNG
    as21 = angle(s21);  % Angle of S21
    das21 = [0; diff(as21)];    % Diff of angle of S21
    [~, pks] = findpeaks(max(abs(das21),6)); 
    val = sign(das21(pks));
    m3 = m2*0;  m3(pks) = val;  m3 = cumsum(m3);
    isDPS = 1-2*(as21 - 2*pi*m3 > 0).*(X < 1);
    m1 = m1.*isDPS; m2 = m2.*isDPS;

    % Recalculate Log Branch (and DNG) Corrected C-parameters
    [zr, n, pro, X, Yn, Yd] = CParam(f, s11, s21, d, m1, m2);
    
    % Correct sign of zr for Q2
    zr = zr.*(1-2*( (X>1) & (imag(s11) > 0) ));

    u = n.*zr;
    e = n./zr;
    
end