%
%  This routine determines the controller for the quadratic optimal system
%
%  Inputs: Gp = plant transfer function
%           q = weight
%  Outputs: Go = desired closed loop transfer function  or [] for errors
%
function Go = solve_quadratic(Gp,q)
%
  [num_Gp,den_Gp] = tfdata(Gp,'v');
%  
%  first compute D(s)*D(-s)
%
  Dp = den_Gp;
  Nden = length(Dp);
  entry = 1;
  scale = [entry];
  for i = 2:Nden
      entry = -entry;
      scale = [scale entry];
  end;
  Dm = Dp.*fliplr(scale);
  DD = conv(Dp,Dm);
%
%  Now compute N(s)*N(-s)
%
  Np = num_Gp;
  Nnum = length(Np);
  entry = 1;
  scale = [entry];
  for i = 2:Nnum
      entry = -entry;
      scale = [scale entry];
  end;
  Nm = Np.*fliplr(scale);
  NN = conv(Np,Nm);
%
%  Put it all together
%
  Q = DD + q*NN;
%  
%  Now find the roots of this
%
  r = roots(Q);
%
%  Now we need to find the stable roots (negative real parts)
%
  stable = [];
  for i = 1:length(r)
      if (real(r(i)) < 0 ) 
          stable = [stable r(i)];
      end;
  end;
%
%  Now construct D0, the denominator of the transfer function
%
  D0 = poly(stable);
%
%  Now put it all together
%
  num_Go = (q*num_Gp(end)/D0(end))*num_Gp;
  den_Go = D0;
  Go = tf(num_Go,den_Go);
%
  return;