%
  disp(' ')
  disp(' -------------Open Loop Linear Analysis:-------------')
  disp(' ')
%
% Form the mass matrices:
  MSway  = [  mass-Yvdot,     mass*xG-Yrdot;
              mass*xG-Nvdot,        Izz-Nrdot];
  MSurge =    mass-Xudot;
%
% Form the AA matrices:
  YrudFixed = .5*rho*ControlArea*vel*(Cd0+dCl);
  AASway  = [  Yv*vel-YrudFixed,          Yr*vel-mass*vel-xR*YrudFixed;
               Nv*vel-xR*YrudFixed,   Nr*vel-mass*xG*vel-xR^2*YrudFixed];
%
% AASwayIdeal has no rudder, just a pure turning moment
  AASwayIdeal = [  Yv*vel,      Yr*vel-mass*vel;
                   Nv*vel,   Nr*vel-mass*xG*vel];
  AASurge = (2*Xuu-rho*ControlArea*Cd0)*vel;
%
% Form the control matrices.  The first term in Ydelta is from the rudder,
% and the second is the contribution from the vectored thrust.
  Ydelta = .5*rho*ControlArea*vel^2*dCl;
  BBSway  = [    Ydelta + Tp;
              xR*Ydelta + xP*Tp];
  BBSurge = 1;
%
% The ideal system has no rudder, just a pure turning moment of xR*Yideal Nm.
  Yideal  = 75;
  BBSwayIdeal  = [  Yideal;
                  xR*Yideal];
%
% If Ideal System;
%  AASway = AASwayIdeal;
%  BBSway = BBSwayIdeal;
%
% First, compute the steady-state values.  Accelerations are zero.
  ySwaySS  = -AASway\BBSway*deltaR;
  betaCGSS = ySwaySS(1,1)/vel;
  Rturn    = abs(vel/ySwaySS(2,1));
  Ncontrol = ( xR*Ydelta + xP*Tp )*deltaR;
  NcontrolP= xPcp*Ydelta*deltaR;
  NrudderSS= Ncontrol - xR*YrudFixed*(ySwaySS(1,1)+ySwaySS(2,1)*xR);
  NpivotSS = NcontrolP-xPcp*YrudFixed*(ySwaySS(1,1)+ySwaySS(2,1)*xR);
  disp(' Steady-State Values of the OPEN LOOP Linear Model:')
  disp([' The rudder angle is ', num2str(deltaR*r2d), ' Deg.'])
  disp([' The CG drift angle is ', num2str(betaCGSS*r2d),' Deg.'])
  disp([' The turn rate is ', num2str(ySwaySS(2,1)*r2d),' Deg/Sec'])
  disp([' The turn diameter is ', num2str(2*Rturn),' Meters',...
        ' or ', num2str(2*Rturn/len), ' vehicle lengths.'])
  disp([' The control-torque component about the CG is ',...
        num2str(Ncontrol),' N-m'])
  disp([' The control-torque component about the pivot point is ',...
        num2str(NcontrolP),' N-m'])
  disp([' The total steady-state torque about the CG is ',...
        num2str(NrudderSS),' N-m'])
  disp([' The total steady-state torque about the pivot point is ',...
        num2str(NpivotSS),' N-m'])
%
% Form state-space representations:
  ASway = MSway\AASway;
  BSway = MSway\BBSway;
  CSway = eye(2);
  DSway = [0;0];
  SwaySys  = ss(ASway, BSway, CSway, DSway);
  SurgeSys = ss(MSurge\AASurge,  MSurge\BBSurge, 1, 0);
%
% Check the eigenvalues:
  disp(' The eigenvalues are: ')
  evals = eig(ssdata(SwaySys))
  disp(' The natural frequency, damping ratio, and time constant are:')
  disp ' f_n , damping ratio, time constant '
  [abs(evals)/2/pi -real(evals)./abs(evals) -1./real(evals)]
%
% Simulate a turn;
  [ySway, t] = step(SwaySys,tsim);
%
% Scale for rudder deflection:
  ySway = ySway*deltaR;
%
% OK, lets put this baby in a loop and see where the eigenvalues go as
% a function of control surface area:
%
  for i=1:Npts,
    Area(i) = ControlArea*i/Npts*2;
    YrudFixed1  = .5*rho*Area(i)*vel*(Cd0+dCl);
    AASway1  = [  Yv*vel-YrudFixed1,          Yr*vel-mass*vel-xR*YrudFixed1;
                  Nv*vel-xR*YrudFixed1,   Nr*vel-mass*xG*vel-xR^2*YrudFixed1];
    poles(i,:) = eig(MSway\AASway1)';
%
%   Form the control matrices.  The first term in Ydelta is from the rudder,
%   and the second is the contribution from the vectored thrust.
%
    Ydelta1 = .5*rho*ControlArea*vel^2*dCl;
    BBSway1 = [    Ydelta1 + Tp;
                   xR*Ydelta1 + xP*Tp];
%
    ySwaySS1     = -AASway1\BBSway1;
    betaCGSS1    = ySwaySS1(1,1)/vel;
    Rturn1       = abs(vel/ySwaySS1(2,1));
    Ncontrol1    = ( xR*Ydelta1 + xP*Tp );
    NrudderSS1(i)= Ncontrol1 - xR*YrudFixed1*(ySwaySS1(1,1)+ySwaySS1(2,1)*xR);
  end
  Area = Area';
% NrudderSS1 is the deltaR to Nm about the CG ss value.
  disp(' ')
  disp(' -------------Closed Loop Linear Analysis:-------------')
  disp(' ')

%
% First, form the full open loop Sway system including the linearized 
% kinematics:
  SwaySysOl = ss([ASway    0*eye(2);
                   eye(2) [0   vel;
                           0    0]] ,...
                 [BSway'   0    0 ]',...
                      eye(4),...
                 [0  0     0    0 ]');
  disp(' The full order (w/ kinematics) open-loop sway system e-values are:')
  evals = eig(ssdata(SwaySysOl))
%