% RecipCal.m
% Matlab/Octave script to calculate solution to reciprocity calibration with syms

pkg load symbolic % necessary when running in Octave

syms s T Trev Tv Ze To Zm Cm ZinOpen ZinShort ToSq 

T=[[1 Ze];[0 1]]*[[0 To]; [1/To 0]]*[[1 Zm]; [0 1]]; % SRR brackets around rows to avoid a bug in Python
DetT=det(T);
Trev= (1/DetT)*[[T(2,2) T(1,2)]; [T(2,1) T(1,1)]]; % T with ports reversed
Tv=[[1 0]; [s*Cm 1]]; % T for volume between speakers
M=simplify((T*Tv*Trev)); 
M = subs(M,To^2,ToSq); % ToSq= To^2

shouldBeZero = simplify(M(1,1)-M(2,2)); % its is 0 (confirms A=D, reversible) 

ZinOpen_Alg = simplify(M(1,1)./M(2,1)); %A/C
ZinShort_Alg = simplify(M(1,2)./M(2,2)); %B/D

ToSq_AlgOpen = simplify(solve(ZinOpen_Alg==ZinOpen,ToSq)); % single valued
ToSq_AlgShort = simplify(solve(ZinShort_Alg==ZinShort,ToSq)); % multivalued

Zm_AlgOpen = simplify(solve(ZinOpen_Alg==ZinOpen,Zm)); % multivalued
Zm_AlgShort = simplify(solve(ZinShort_Alg==ZinShort,Zm)); % multivalued

% Substitute ToSq expression (since it is single valued), then solve for Zm
ZinShort_Alg_Zm = simplify(subs(ZinShort_Alg,ToSq,ToSq_AlgOpen));
Zm_sol = solve(ZinShort_Alg_Zm == ZinShort,Zm);

% Sarah: Students should use 
% Zm_sol: Solution for Zm in terms of known parameters 
% ToSq_AlgOpen: Solution for To^2 in terms of Zm and known parameters

% construct Matlab vector equations, allow for '.' vector operations
Zm_eq1 = char(Zm_sol(1,:));
Zm_eq1= strrep(Zm_eq1,'*','.*'); Zm_eq1= strrep(Zm_eq1,'/','./');  
Zm_eq2 = char(Zm_sol(2,:));
Zm_eq2= strrep(Zm_eq2,'*','.*'); Zm_eq2= strrep(Zm_eq2,'/','./');  
ToSq_eq = char(ToSq_AlgOpen);
ToSq_eq= strrep(ToSq_eq,'*','.*'); Zm_eq= strrep(ToSq_eq,'/','./');  

% Results:
% (Choose the Zm that is "positive real)
% Zm_eq1 = (-ZinOpen.*(ZinOpen - ZinShort) + sqrt(ZinOpen.*(ZinOpen - ZinShort)).*(-Ze + ZinOpen))./(Cm.*ZinOpen.*s.*(ZinOpen - ZinShort));
% Zm_eq2 = (-ZinOpen.*(ZinOpen - ZinShort) - sqrt(ZinOpen.*(ZinOpen - ZinShort)).*(-Ze + ZinOpen))./(Cm.*ZinOpen.*s.*(ZinOpen - ZinShort))
% ToSq_eq = -Zm.*(Cm.*Ze.*Zm.*s - Cm.*ZinOpen.*Zm.*s + 2.*Ze - 2.*ZinOpen)/(Cm.*Zm.*s + 1);


