%CParam.m
%[zr, n, pro, X, Yn, Yd] = CParam(f, s11, s21, d, m1, m2)
% Written by Jacob Lewis April 2020.
% This code is to derive the Constitutive Parameters from the Scattering
% Parametres and some other data.
%
% Inputs:
%   f   Frequency
%   s11 Reflection
%   s21 Transmission
%   d   Thickness
%   m1  Sign branch
%   m2  Log branch
% Return:
%   e   Permittivity?
%   u   Permeability?
%   n   Index of refraction
%   zr  Relative Intrinsic Impedance
%   egd Wave propagation < e^(j \gamma d) >
%   X   Variable used in calculating egd <egd = X ± sqrt(X^2 - 1)>

function [zr, n, pro, X, Yn, Yd] = CParam(f, s11, s21, d, m1, m2)
    % Define constant and w
    c = 299792458;
    w = 2*pi*f;
    
    % Solve Stuff
    X = (1 + s21.^2 - s11.^2)./s21/2;
    % If there will be loss of information due to the largness of X^2
    if any(X > 1e7)
        X = sym((1 + s21.^2 - s11.^2)./s21/2); % Convert to symbol
        % Use Symbolic Math
        pro = double(real(X) + m1.*sqrt(real(X).^2 - 1));
        X = double(X);  % Convert back to double
    else
        pro = real(X) + m1.*sqrt(real(X).^2 - 1);
    end
    n = (-1j*log(pro) + 2*pi*m2)*c./w/d;
    
    % Two different calculations for if Zr is primarily real or imaginary
    Yn = s21.^2 - (1 + s11).^2;
    Yd = (s21.^2 - (1 - s11).^2);
    zr = sqrt(real(Yn./Yd));
end
