% To use S parameters in Simulink, one must convert the frequency domain
% data into a model that is expressed as differential equations. 
% The differential equations are expressed as Laplace Transforms  typically
% using the variable "s" (as in s domain) which can get confusing since we are also using
% the variable S to represent Scattering Parameters. 
% 
% To get to the s-domain the invfreqs function will be used. Transistors
% have rather low order models, so this simple function works reasonably
% well. 
 
clc
clear
ITER = 40;    % set some stopping parameters for invfreqs
TOL  = 0.01;

warning('off') % makes us feel better ... :-) 
[sinp,f] = s_extract('mrf904_10_10a.rfa',10,500,10);   % load the S parameter data 
% cycle through s11, s12, s21, s22
  for m=1:2
     for n=1:2
      h = sinp(m,n).s;
      w = 2*pi*f;       % note that f is in MHz ... 
      nb  = 4;
      na  = 4;
      run = 1;
      if m==n
         weight=[];
     else
         weight=1./(w);
     end
      [m,n]
      while run
         lastwarn('');
         [s904(m,n).num, s904(m,n).den] = invfreqs(h,w,nb,na,weight,ITER,TOL);     
         wmsg = lastwarn;
         if strcmp('',lastwarn) | (na==2) | (nb==2)
             run = 0;
         else
             nb=max(nb-1,2);
             na=max(na-1,2);
         end;
      end;
      % plot the poles and zeros of the s-plane models
      hf = figure('name',sprintf('s %i %i',m,n));
      plot(roots(s904(m,n).den)+j*eps,'x');      % need the eps to force it to plot real vs imag when there is no imag ! 
      hold on
      plot(roots(s904(m,n).num)+j*eps,'o');
      
      % Compare the measured magnitude and phase with the mag and phase
      % computed from the s-domain models. Although agreement is not
      % perfect, it is quite reasonable. 
      hf = figure('name',sprintf('s %i %i',m,n));
      freqs(s904(m,n).num,s904(m,n).den,w)
      hax=get(hf,'children');
      line('parent',hax(1),'xdata',w,'ydata',abs(h),'color',[1 0 0])
      line('parent',hax(2),'xdata',w,'ydata',(180/pi)*atan2(imag(h),real(h)),'color',[1 0 0])
   end;
 end;
 
 warning('on')             % restore the warning state
 save mrf904_smodel  s904  % dave the s-domain results for future use in Simulink