% MassCal.m
% Matlab/Octave script to calculate Hunt Parameters with syms

pkg load symbolic % necessary when running in Octave

% Calculate Hunt parameters using 3 input impedance measurements (Zin), 3 load impedance models (ZL)
syms Zm Ze ToSq ZinL1 ZinL2 ZinL3 ZL1 ZL2 ZL3 

% 3 equations in 3 unknowns (these assume an anti-reciprocal system)
ZinL1_Alg = Ze + ToSq/(ZL1 + Zm);
ZinL2_Alg = Ze + ToSq/(ZL2 + Zm);
ZinL3_Alg = Ze + ToSq/(ZL3 + Zm);

% construct 2 equations in Zm & ToSq (eliminate Ze)
EQ1 = simplify(ZinL1_Alg -ZinL2_Alg);
EQ2 = simplify(ZinL2_Alg - ZinL3_Alg);

ToSq_Zm = simplify(solve(EQ1 == (ZinL1-ZinL2),ToSq)); % Find To^2 in terms of Zm_sol
EQ2_Zm = simplify(subs(EQ2,ToSq,ToSq_Zm)); % Sub this expression for To^2 into the other equation
Zm_sol = (solve(EQ2_Zm==(ZinL2-ZinL3),Zm)); % Solve for Zm
ToSq_sol = simplify(subs(ToSq_Zm,Zm,Zm_sol)); % Solve for To^2

ZinL1_Alg2 = subs(ZinL1_Alg,Zm,Zm_sol); % Sub Zm into EQ 1
ZinL1_Alg3 = subs(ZinL1_Alg2,ToSq,ToSq_sol); % Sub To^2 into EQ 1
Ze_sol = simplify(solve(ZinL1_Alg3==ZinL1,Ze)); % Solve for Ze

% Simplify To^2 expression by factoring numerator and denominator
[nToSq_sol,dToSq_sol] = numden(ToSq_sol);
ToSq_sol = factor(nToSq_sol)/factor(dToSq_sol);

% construct Matlab vector equations, allow for '.' vector operations
Zm_eq = char(Zm_sol);
Zm_eq= strrep(Zm_eq,'*','.*'); Zm_eq= strrep(Zm_eq,'/','./');  
Ze_eq = char(Ze_sol);
Ze_eq= strrep(Ze_eq,'*','.*'); Ze_eq= strrep(Ze_eq,'/','./'); 
ToSq_eq = char(ToSq_sol);
ToSq_eq= strrep(ToSq_eq,'**2','.^2'); % there is a bug: ** = ^
ToSq_eq= strrep(ToSq_eq,'*','.*'); ToSq_eq= strrep(ToSq_eq,'/','./'); 

% Results:
% Zm_eq = (ZL1.*ZL2.*ZinL1 - ZL1.*ZL2.*ZinL2 - ZL1.*ZL3.*ZinL1 + ZL1.*ZL3.*ZinL3 + ZL2.*ZL3.*ZinL2 - ZL2.*ZL3.*ZinL3)./(ZL1.*ZinL2 - ZL1.*ZinL3 - ZL2.*ZinL1 + ZL2.*ZinL3 + ZL3.*ZinL1 - ZL3.*ZinL2)
% Ze_eq = (ZL1.*ZinL1.*ZinL2 - ZL1.*ZinL1.*ZinL3 - ZL2.*ZinL1.*ZinL2 + ZL2.*ZinL2.*ZinL3 + ZL3.*ZinL1.*ZinL3 - ZL3.*ZinL2.*ZinL3)./(ZL1.*ZinL2 - ZL1.*ZinL3 - ZL2.*ZinL1 + ZL2.*ZinL3 + ZL3.*ZinL1 - ZL3.*ZinL2)
% ToSq_eq = -(ZL1 - ZL2).*(ZL1 - ZL3).*(ZL2 - ZL3).*(ZinL1 - ZinL2).*(ZinL1 - ZinL3).*(ZinL2 - ZinL3)./(ZL1.*ZinL2 - ZL1.*ZinL3 - ZL2.*ZinL1 + ZL2.*ZinL3 + ZL3.*ZinL1 - ZL3.*ZinL2).^2