% file: PRtrans.m
%     for the transformation of a minimum degree two PR functions state 
%     equations to a passive coupling circuit.
%  RWN  11/22/06  uses the CARE program to solve the Riccati equation
%  for y(s)=(s^2+as+b)/(s^2+cs+d)
%------------------- set variables for passing to functions -------
global A1 B1 Do Q Yc Yosym Ycskw
%-------------------
%--------------- input data
a=input('numerator coefficient of s a ')
%     
b=input('numerator constant coefficient b ')
%
c=input('denominator coefficient of s c ')
%
d=input('denominator constant coefficient d ') 
%
%-------------- form state variable matrices
A=[-c 1;-d 0]
B=[a-c;b-d]
C=[1 0]
D=[1]
%
%-------------- form initial coupling Yc -----
Yci=[1 C;-B -A]
Ycisym=(Yci+Yci')/2
%
%------------- form matrices for Riccati equation -
Do=sqrtm((D+D'))
A1=A-B*[Do^(-2)]*C
B1=B*[Do^(-1)]
Q=C'*[Do^(-2)]*C
R=[-eye(1)]
%
%-------------solve the Riccati equation -----------
%  A1'X+XA1-XB1RB1'X+Q=0	
 [X,L,G]=care(A1,B1,Q,R)
%-----------form PR YC ---------
	 T1=sqrtm(X)	 
	 Yc=[D C*T1^(-1);-T1*B -T1*A*T1^(-1)]
	 Ycsym=(Yc+Yc')/2
         Ycskw=(Yc-Yc')/2
%----------- factor the symmetric part
%	 Yo=chol((Yc+Yc')/2)
%        Yo*Yo'
% --------- due to numerical error chol(.) does not work; use SVD
         [U,S,V]=svd(Ycsym)
% gives Ycsym=U*diagS*V'
%--------- final decomposition for coupling Y symmetric part
          Yosym=U(:,1)*S(1,1)*U(:,1)'
% Yosym is realized by one resistor seen through a transformer of
% turns ratio matrix U(:,1)
% final synthesis is a parallel connection of Ycskew with Yosym 
%  loaded in two capacitor at the last two ports. 

end