%--------------------------------------------------------------------------
%
% PURPOSE: Simulate the steering dynamics of Altex for small angles,
%          in the horizontal plane.
% AUTHOR:  Rob McEwen
% DATE:    99/2/12
%
%--------------------------------------------------------------------------
%
% Assume that: 1) The vehicle is in straight and level flight, so this
%                 analysis is in the horizontal plane.
%              2) All angles remain less than 6 degrees.
%              3) The ring and struts have the lift and drag of 4 
%                 Odyssey fins.
%
% NOTES: 1) The effect of the tailcone body deflection on the hydrodynamic 
%           derivatives is not modeled here.
%        2) The signs of Nv and Yr are flipped from Knut's eqns due to
%           a suspected reference frame flip in Landweber et. al.
%
% Units and Fundamental Constants:
  in2m = .0254;             % Inches to meters.
  d2r  = pi/180;            % Degrees to radians
  r2d  = 180/pi;            % Radians to degrees
  rho  = 1028;              % Density of sea water at 0 Deg C; kg/m^3
  s2kg = 14.59390;          % Slugs to kg
  ft2m = .3048;             % Feet to meters
  sft2kgm2 = 1.355818;      % Slug-ft^2 to kg-m^2
  lbf2n    = 4.448;         % Pounds force to Newtons
  kph2mps  = 1852/3600;     % Knots per hour to meters per sec
  ftlbf2nm = 1.3558179;     % Foot-pounds force to Newton-meters.
  oi2nm    = .0070615516;   % Ounce-inches to Newton-meters.
%
  caseName = '-1Deg';        % Plot name suffix
  caseName = 'Tp52Dr10';     % Plot name suffix
%
% Altex/Dorado Specific Constants:
  vel = 3.0*kph2mps;        % speed in m/s
  Tp  = 52;                 % Thrust in Newtons
%  Tp = 0;
  mass= 795;                % Mass; kg
  len = 4.19;               % Length; meters
  dia = 21*in2m;            % Diameter; meters
  xG  = 0;                  % Longitudinal CG offset, meters
  xPcp= -5*in2m;            % Distance from pivot to ring CP in tailcone coords
  aR  = 3.25*in2m;          % Dist from c.l. to rudder actuator base
%
% Altex/Dorado constants taken from the Odyssey sim:
% NOTE: Assume that for Yaw/Pitch control purposes, the ring has 
% the properties of an Odyssey fin, with 4 times the surface area.  For 
% the surge equation, use the actual (approximate, of course) area.
% 
  StrutArea  = 0.0299;      % Odyssey fin area; m^2
  ControlArea = StrutArea*4; % Equivalent rudder/elevator area (m^2). 
%
  RingArea = 24*in2m*pi*10*in2m + StrutArea*4;
                         % Ring + 4 struts;m^2. For drag calc in the Surge Eqn.
  Ar  = 4.26;            % Aspect ratio
  Cdc = 0.723;           % Crossflow Cd
  Cd0 = 0.0065;          % Cd from NACA 0015 section; PNA p. 503
  Ef  = 0.90;            % "Oswald" efficiency factor; PNA p. 503
  dCl = 3.392;           % Lift coefficient slope 
  Xu  = 10.0;            % Body drag; N/(m/s) at "vel" m/s
%
% Computed Constants:
  xR  = -len*.45;           % Distance from CG to rudder; m
  xP  = xR - xPcp;          % Distance from CG to pivot point.
  Izz = mass/5*((len/2)^2+(dia/2)^2); %Ellipsoid; Greenwood p501. kg-m^2
  Ixx = 2*mass/5*(dia/2)^2; % Greenwood p501. Assume homogenous dist. kg-m^2
  Nr  = -.045*mass*len;     % Deriv of N wrto r; kg-m (*m/sec^2=Nm), L&J
  Nv  = -.795*mass;         % Deriv of N wrto v; kg (*m^2/sec^2=Nm), L&J
  Yr  = .064*mass;          % Deriv of Y wrto r; kg (*m/sec^2=N), L&J
  Yv  = -.5*mass/len;       % Deriv of Y wrto v; kg/m   (*m^2/sec^2=N), L&J
  Xudot = -.1*mass;         % Deriv of X wrto u-dot (Added Mass); kg; WAG
  Yvdot = -.8*mass;         % Deriv of Y wrto v-dot; (Added Mass) kg; WAG
  Nrdot = -.8*Izz;          % Deriv of N wrto r-dot; (Added Iner) kg-m^2; WAG
  Yrdot = 0;                % Deriv of Y wrto r-dot; (Coupled) kg-m; WAG
  Nvdot = 0;                % Deriv of N wrto v-dot; (Coupled) kg-m; WAG
%
% Nonlinear terms:
  Yvabsr = 0;               % Knut
  Yvabsv = -2.94*mass/len;  % Deriv of Y wrto v|v|; kg/m; Knut
  Nrabsv = 0;               % Knut
  Nvabsv = .344*mass;       % Deriv of N wrto v|v|; kg; Knut
  Nrabsr = -.0705*mass*len^2;% Deriv of N wrto r|r|; kg-m^2; Knut
%
% Sim Parameters:
  tf   = 60;                % Stop time, sec
  Npts = 500;               % Number of points;
  deltaR = 1*d2r;
  deltaR = 10*d2r;
  tdeltaR= 2;               % Rudder slew time
  tdur= 20;                 % Rudder slew time
%
% Flags:
%
% VehFlg = 0 => Altex;
% VehFlg = 1 => Odyssey;
% VehFlg = 2 => DSRV;
  if exist('hardflag')==0, hardflag=0, end
  if exist('VehFlg')==0, VehFlg=0, end
  LinFlg = 0;                           % 1=> Plot linear as well; 0=> sim only
%
% Start screen log:
  !rm run.txt
  diary(['run',caseName,'.txt'])
  disp(' ')
  if VehFlg == 0,  disp(' Altex Parameters:'), end
%
%----------START  ODYSSEY---------------
%
% Lets try it with Odyssey numbers:
  if VehFlg == 1,
    disp(' Using Odyssey Parameters:')
    tf   = 20;                % Stop time, sec
%   vel = 1852*3.0/3600;      % speed in m/s ????
    mass= 339;                % Mass; kg
    len = 4.19/2;             % Length; meters  ????
    ControlArea = StrutArea*2; % Equivalent rudder/elevator area (m^2). 
    xR  = -1.069;             % Distance from CG to rudder; m
    Izz = 71.95;              % kg-m^2
    Ixx = 10.59;              % kg-m^2
    Nr  = -74.9;              % Deriv of N wrto r; kg-m^2/sec (*1/sec=Nm), L&J
    Nv  = -167.3;             % kg-m/sec
    Yr  = 112.6;              % kg-m/sec
    Yv  = -117.5;             % kg/sec
    Yvdot = -313.6;           % kg
    Nrdot = -50.7;            % kg-m^2
    Yvabsr = 0;               
    Yvabsv = -4930;           % Deriv of Y wrto v|v|; kg/m; sim.c
    Nrabsv = 0;               
    Nvabsv = 0;               % Deriv of N wrto v|v|; kg; sim.c
    Nrabsr = -150;            % Deriv of N wrto r|r|; kg-m^2; sim.c
  end
%
%----------END ODYSSEY---------------
%----------START  DSRV---------------
%
% Lets try it with DSRV numbers.  These are all from pp. 36,37 & 46.
  if VehFlg == 2
    disp(' Using DSRV Parameters:')
    tf  = 100;                % Sim stop time, sec
%   vel = 1852*3.0/3600;      % speed in m/s
    mass= 4363*s2kg;          % Mass; kg
    len = 49.33*ft2m;         % Length; meters 
    ControlArea = 8*1.83*(ft2m)^2; % Equivalent rudder/elevator area (m^2). 
    xR  = -(46.89-23.4)*ft2m;% Distance from CG to rudder; m
    Izz = 475868*sft2kgm2;    % kg-m^2
    Ixx = 38244*sft2kgm2;     % kg-m^2
    Nr  = -49200*sft2kgm2;    % Deriv of N wrto r; kg-m^2/sec (*1/sec=Nm)
    Nv  = -3050*s2kg*ft2m;    % kg-m/sec
    Yr  = 792*s2kg*ft2m;      % kg-m/sec
    Yv  = -136.2*s2kg;        % kg/sec
    Yvdot = -4030*s2kg;       % kg
    Nrdot = -351000*sft2kgm2; % kg-m^2
    Yvabsr = 0;               % 
    Yvabsv = -2.38;           % Deriv of Y wrto v|v|; kg/m;
    Nrabsv = 0;               % 
    Nvabsv = 0;               % Deriv of N wrto v|v|; kg;
    Nrabsr = 0;               % Deriv of N wrto r|r|; kg-m^2;
    Cd0 = 0.05;               % 
    dCl = 1.35/(15*pi/180);   % Lift coefficient slope 
    Cdc =    0;               % Crossflow Cd
    Ef  = 0.90;               % "Oswald" efficiency factor; PNA p. 503
    Ar  = (20*d2r*dCl)^2/.35/pi/Ef;   % Aspect ratio
  end
%
%----------END DSRV---------------
%
% Print parameters to the screen and diary file.
  disp(' ')
  document
  disp(' ')
  disp(' -------------Linear Analysis:-------------')
  disp(' ')
%
% Set up time vector:
  tsim = [0:Npts]'*tf/Npts;
%
% 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 = Xu-rho*RingArea*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 = Tp - .5*rho*RingArea*vel^2*Cd0 - Xu*vel;
%
% 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 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 representation:
  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.





%
% -------------Simulink--------------
%
  disp(' ')
  disp(' -------------Simulation of Nonlinear EOM:-------------')
  disp(' ')
%
% Build parameter vector:
  MSwayInv = inv(MSway);
%
% DEBUG: Zero out the nonlinear terms:
%  Cdc = 0;
%  Yvabsv = 0;
%  Nvabsv = 0;
%  Nrabsr = 0;
%
  CDynVec = [ MSwayInv(1,:)';
              MSwayInv(2,:)';
              mass;
              vel;
              xR;
              xP;
              aR;
              dCl;             % Coeff of lift; linear coeff
              Cdc/Ar;          % Coeff of lift; second order coeff
              Cd0;             % Coeff of drag; constant term
              1/pi/Ar/Ef;      % Coeff of drag; Cl^2 coeff
              .5*rho*ControlArea; % Coeff for L and D eqns
              Yv;
              Yr;
              Yvabsr;
              Yvabsv;
              xG;
              Nv;
              Nr;
              Nrabsv;
              Nvabsv;
              Nrabsr;
            ]';
%
% Run it!
  [tout] = sim('AltexSim',tsim,simset('maxstep',tsim(2)-tsim(1)));
%
% -------------End Simulink--------------
%
%
% Indirect addressing for the columns of yout.
  lbls = { 'vIx'
           'rIx'
           'betaIx'
           'betaRIx'
           'AofAIx'
           'TrqNRIx'
           'TrqPIx'
           'xPosIx'
           'yPosIx'
           'YawIx'
           'FYRIx'
           'FYTIx'
           'FRAIx'
           'deltaRIx'
         };
  sz = size(lbls);
  i = 1;
  while( i <= sz(1) ),
    eval([ char(lbls(i))  ' = '  num2str(i) ';']);
    i = i + 1;
  end
%
%
% ------------ Nonlinear Steady-State Analysis --------------
%
% Now compute some steady-state values from the nonlinear steering
% equations from Knut's derivation.  See my notes, p. 12, 99/3/2.
% I've assumed that rCmd < 0, so this won't work for positive values 
% without a tweak.
%
% WARNING: This will match the sim SS values ONLY IF Tp = 0.
%
% Also, the deltaR from above does not apply.  The deltaR here is different.
%
  rCmd = yout(Npts+1,rIx);
  Ap = (xR*Yvabsv - Nvabsv)/mass/vel^2;
  Bp = (xR*Yv - Nv)/mass/vel;
  Cp = ((xR*Yr - Nr)/mass - xR)*rCmd/vel + Nrabsr/mass*(rCmd/vel)^2;
  vSS = (-Bp + sqrt(Bp^2-4*Ap*Cp))/2/Ap;  %Positive root of a quadratic
  FYRSS = mass*vel*rCmd - Yr*vel*rCmd - Yv*vel*vSS - Yvabsv*vSS^2;
  turnDia = 2*vel/abs(rCmd);
  disp(' ------- Nonlinear Steady-State Analysis: ----------')
  disp(' ')
  disp(sprintf(' %-15f: Computed vSS, Meters/Sec', vSS ));
  disp(sprintf(' %-15f: Simulated vSS, Meters/Sec', yout(Npts+1,vIx)));
  disp(sprintf(...
       ' %-15f: Tp. MUST BE ZERO, or else this analysis does not apply.', Tp));
  disp(sprintf(' %-15f: Yaw rate, Deg/Sec.', rCmd/d2r));
  disp(sprintf(' %-15f: Turn Diameter, Meters.', turnDia));
  disp(sprintf(' %-15f: Beta, Deg/Sec.', asin(vSS/vel)/d2r));
  disp(sprintf(' %-15f: BetaR, Deg/Sec.', asin((vSS + xR*rCmd)/vel)/d2r));
  disp(sprintf(' %-15f: Rudder force resolved to body y axis, Newtons.',FYRSS));
  disp(sprintf(' %-15f: Simulated Rudder force, Newtons.',yout(Npts+1,FYRIx)));
%
  diary off    
  clear functions mex
  str = [', Case:',caseName,' RSM ',myDate];
  PaperPosn = [.5 .5 3.5 2.5];

  figure(1)
  if hardflag ==3, set(gcf,'PaperPosition',PaperPosn), end
  clf
  subplot(2,1,1)
  if LinFlg, 
    plot(t, [ySway(:,1) yout(:,vIx)]);
  else
    plot(t, [yout(:,vIx)]);
  end
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  ylabel('Sway Rate')
  title(['Sway Rate, Meters per Second',str]);
  grid on
  if LinFlg, legend('Linear','NonLin'), end
  subplot(2,1,2)
  if LinFlg, 
    plot(t, [ySway(:,2) yout(:,rIx)]/d2r);
  else
    plot(t, [yout(:,rIx)]/d2r);
  end
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  ylabel('Yaw Rate, Degrees/Second')
  title(['Yaw Rate, Degrees/Second',str]);
  xlabel('Seconds');
  if LinFlg, legend('Linear','NonLin'), end
  grid on
%
% Compute drift angles at CG and stern.  See notes, Eqn.'s 8 & 9.
  figure(2)
  if hardflag ==3, set(gcf,'PaperPosition',PaperPosn), end
  clf
  beta  = ySway(:,1)/vel;
  betaR = (ySway(:,1)+xR*ySway(:,2))/vel;
  AofA  = (deltaR-betaR);
  beta_NL  = yout(:,1)/vel;
  betaR_NL = (yout(:,1)+xR*yout(:,2))/vel;
  AofA_NL  = (deltaR-betaR);
  if LinFlg,
    plot(t, [beta yout(:,betaIx), betaR, yout(:,betaRIx),...
             AofA, yout(:,AofAIx)]*r2d)
    legend('Lin CG Drift', 'Nonlin CG Drift', 'Lin Stern Drift', ...
           'Nonlin Stern Drift', 'Lin Rudder Attack', ...
           'Nonlin Rudder Attack',1)
  else
    plot(t, ...
    [yout(:,betaIx), yout(:,betaRIx), yout(:,AofAIx) yout(:,deltaRIx)]*r2d)
    legend('Drift Angle at CG', 'Drift Angle at Stern',...
           'Angle of Attack on the Rudder','deltaR',1)
  end
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  title(['Drift Angle at CG, Drift Angle at Rudder, Angle of Attack',str])
  ylabel('Degrees')
  xlabel('Seconds');
  grid on
%
% Compute turning moment due to rudder.  See notes, Eqn. 32.  Also
% compute the torque about the tailcone pivot point.
%
  figure(3)
  if hardflag ==3, set(gcf,'PaperPosition',PaperPosn), end
  clf
  FYR      = Ydelta*deltaR - YrudFixed*(ySway(:,1)+ySway(:,2)*xR);
  Nrudder  = xR*FYR;
  Nturn    = Nrudder + xP*Tp*deltaR;
  Npivot   = xPcp*(Ydelta*deltaR - YrudFixed*(ySway(:,1)+ySway(:,2)*xR));
  if LinFlg,
    plot(t, [Nturn yout(:,TrqNRIx) Npivot yout(:,TrqPIx)]);
    legend('Lin Vehicle CG','Nonlin Veh CG','Lin Pivot','Nonlin Pivot',0)
  else
    plot(t, [yout(:,TrqNRIx) yout(:,TrqPIx)]);
    legend('Torque About Veh CG','Torque About Pivot Pt',0)
  end
%  if LinFlg,
%    plotyy(t, [Nturn yout(:,TrqNRIx) Npivot yout(:,TrqPIx)],...
%           t, [Nturn yout(:,TrqNRIx) Npivot yout(:,TrqPIx)]/ftlbf2nm);
%    legend('Lin Vehicle CG','Nonlin Veh CG','Lin Pivot','Nonlin Pivot',0)
%  else
%    plotyy(t, [yout(:,TrqNRIx) yout(:,TrqPIx)],...
%           t, [yout(:,TrqNRIx) yout(:,TrqPIx)]/ftlbf2nm);
%    legend('Torque About Veh CG','Torque About Pivot Pt',0)
%  end
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  title(['Moments about CG and Tailcone Pivot Point',str])
  ylabel('Newton-Meters')
  xlabel('Seconds');
  grid on
%
  figure(4)
  if hardflag ==3, set(gcf,'PaperPosition',PaperPosn), end
  clf
  if LinFlg,
    plot(t, [FYR yout(:,FYRIx)]);
    legend('Lin Rud Frc','Nonlin Rud Frc',0)
  else
    plot(t, [yout(:,[FYRIx FYTIx FRAIx])]);
    legend('Rudder on Body Y','Thrust on Body Y','Actuator Axis',0)
  end
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  title(['Rudder Force Resolved onto the Vehicle y & Actuator Axes',str])
  ylabel('Newtons')
  xlabel('Seconds');
  grid on
%
  figure(5)
  if hardflag ==3, set(gcf,'PaperPosition',PaperPosn), end
  clf
  plot(yout(:,yPosIx), yout(:,xPosIx))
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  title(['Vehicle Position, Meters for a +',num2str(deltaR/d2r),...
         ' Deg Tailcone Angle',str])
  ylabel('x Position in EF Frame, Meters');
  xlabel('y Position in EF Frame, Meters');
  axis equal
  grid on
%
  figure(6)
  if hardflag ==3, set(gcf,'PaperPosition',PaperPosn), end
  clf
  plot(Area(200:500),-1./poles(200:500,:))
  if hardflag == 3, set(gca,'FontUnits','Normalized'), end
  title(' System Time Constants Versus Rudder Area')
  ylabel('Seconds');
  xlabel(['Meters^2',str]);
  grid on
%
  if hardflag==0,
    return
  elseif hardflag==1,
    optstr = ' -depsc';      %Encapsulated Postscript
  elseif hardflag==2,
    optstr = ' -djpeg90';    %JPEG
  elseif hardflag==3,
    optstr = ' -dtiff';      %TIFF
  else
    disp(' Hardflag must be 0, 1, 2 or 3');
    return
  end

  nfig = 6;
  for i = 1:nfig
    eval(['print -f', num2str(i), optstr, ' simfig',num2str(i), caseName])
  end
  hardflag = 0;

  return



%
%  This was an attempt to compute AofASS.  But, I later realized that
%  it is necessary to find the steady-state Drag since AofASS depends on it.
%  So, its a little more complicated than it first seemed.
%
%  Solve for a deltaR, given rCmd = rSS.
%  First,Find the magnitude of the velocity wrto N at the rudder.
   magVelR2 = vel^2 + (vSS + xR*rCmd)^2;
   magVelR  = sqrt(magVelR2);
   betaRSS = asin( (vSS + xR*rCmd)/magVelR );  % Drift angle at the rudder */
   halfRhoS = .5*rho*ControlArea;
%
%  Assume that AofASS is < 0; so take the negative root.   
   rts = roots([Cdc/Ar -dCl FYRSS/halfRhoS/magVelR2]);
   if rts(1) < 0 
     AofASS = rts(1);
   else
     AofASS = rts(2);
   end
   deltaR = AofASS + betaRSS


