% This script performs a "sanity check" on the modeling done in sp2sd.m and
% the simulation done in the Simulink sparam_1.mdl
clear
syms  a1 a2 b1 b2 bs s11 s21 s12 s22  Gs Gl

% input refl coeff: Gs, load refl coeff: Gl
% input source power bs
Gs = 0.5;    % these were the (arbitrary) source and load reflection coefficients 
Gl = 0.5;    % used in sparam_1.mdl 


% solve the flowchart (Simulink block diagram in sparam_1.mdl) using
% Symbolic Math 
eq1 = bs + Gs*b1 - a1;
eq2 = s21*a1 + s22*a2 -b2;
eq3 = Gl*b2 - a2;
eq4 = a2*s12 + a1*s11 -b1;

sol = solve(eq1,eq2,eq3,eq4,a1,a2,b1,b2);
sol.a1 = simple(sol.a1);    
sol.b1 = simple(sol.b1);
sol.b2 = simple(sol.b2);

sol_b2 = sol.b2;
sol_b2

s11_prime = sol.b1/sol.a1
b2_div_a1 = sol.b2/sol.a1

Av = simple(s21*(1+Gl)/((1-s22*Gl)*(1+s11_prime)));    % voltage gain .... 

fmax=750; % same as Simulink model 
[s,f] = s_extract('mrf904_10_10a.rfa',10,fmax,10);  % get measured s parameter data for mrf904

% eval s11 prime to compare with simulink 
for k=1:length(f)
   x(k) = double(vpa(subs(s11_prime,{s11, s12, s21, s22, Gs, Gl},{s(1,1).s(k),s(1,2).s(k),s(2,1).s(k),s(2,2).s(k),0.5,0.5})));
end;
figure
plot(f,20*log10(abs(x)));
set(gca,'xlim',[0,1.28*fmax],'ylim',[-40,0]); grid on;
text('')