function ps3_1(a)
% function ps3_1(a)
% 
% numerical calculations for checking the analytical solution 
% in Problem 3.1, 6.245/Spring 2004

if nargin<1, a=1; end
s=tf('s');                       % a convenient shortcut 
assignin('base','s',s);          % export s to the workspace
assignin('base','a',a);          % export d to the workspace
load_system('ps3_1a');           % load the design model into workspace
[ap,bp,cp,dp]=linmod('ps3_1a');  % extract the LTI model
close_system('ps3_1a');          % close the design model
p=pck(ap,bp,cp,dp);              % re-write plant model in Mu-Tools format
[k,g]=h2syn(p,1,1,2,0);          % design the controller
[ak,bk,ck,dk]=unpck(k);          % get a state space model of the controller
Kn=ss(ak,bk,ck,dk);              % define controller as a standard LTI object
[ag,bg,cg,dg]=unpck(g);          % closed loop system coefficients
G=ss(ag,bg,cg,dg);               % closed loop model
K=[-2*a^2 -2*a];                 % analytical controller gain
x=-sqrt(0.5*(sqrt(a^2+1)+a^2));  % real part of stable eigenvalue
y=-sqrt(0.5*(sqrt(a^2+1)-a^2));  % imaginary part of stable eigenvalue
L=[1/y;-x/y];                    % observer gain
A=[0 1;a^2 0];                   % open loop plant coefficient matrices
B2=[0;1];
C2=[1 0];
aa=A+B2*K+L*C2;                  % controller coefficient matrices
ba=-L;
ca=K;
da=0;
Ka=ss(aa,ba,ca,da);
Pfi=[2*a^3 2*a^2;2*a^2 2*a];
Pse=[-1/y x/y;x/y -x^2/y-y];
B1=[0 0;1 0];
D12=1;
cost=trace(B1'*Pfi*B1)+trace(D12*K*Pse*K'*D12');
disp('Numerical:')
disp(['      optimal gain squared: ' num2str(norm(G)^2)])
tf(Kn)
disp('Analytical:')
disp(['      optimal gain squared: ' num2str(cost)])
tf(Ka)