function [Ze,ToSq,Zm] = calcHuntParams(Zin,ZL)
 
 % function to calculate Hunt parameters 
 
 % Inputs:
 % Zin - Input impedances (3)
 % ZL -  Load impedances (3)
 
 % Outputs:
 % Ze = Z11
 % ToSq = Z12*Z21
 % Zm = Z22
   
 % fix orientation if necessary
 if size(Zin,1)<size(Zin,2), Zin =Zin.'; end
 if size(ZL,1)<size(ZL,2), ZL =ZL.'; end
   ZinL1 = Zin(:,1); ZinL2 = Zin(:,2); ZinL3 = Zin(:,3);
   ZL1 = ZL(:,1); ZL2 = ZL(:,2); ZL3 = ZL(:,3);
 
 Ze = (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);
        
 Zm = (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);

 ToSq = (ZinL1-Ze).*(Zm+ZL1);
 
