%--------------------------------------------------------------------------
%
% 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.
  m2c= 200*12/in2m;         % Meters-to-counts
  c2m= 1/m2c;               % Counts-to-meters
%
  caseName = '-1Deg';       % Plot name suffix
  caseName = 'ClosedLoop';  % Plot name suffix
%
% Altex/Dorado Specific Constants:
  vel = 3.0*kph2mps;        % speed in m/s.  
  Tp  = 52;                 % Thrust in Newtons
  mass= 795;                % Mass; kg
  len = 4.19;               % Length; meters
  dia = 21*in2m;            % Diameter; meters
  xG  = 0;                  % Longitudinal CG offset, meters
  Decouple = 0;             % Decouple the Surge EOM from Sway/Yaw
  Ts = .1;                  % Control System sampling rate, sec.    
  vel_C_N_N_x = 0;          % Vel of the Current (C) wrto an Earth-Fixed 
  vel_C_N_N_y = 1*kph2mps;  %% Newtonian frame (N), coord in the N frame.
  r_No_WP1_N_x = 500;       % Dist from Earth-Fixed origin to first way point
  r_No_WP1_N_y = 500;       % Dist from Earth-Fixed origin to first way point
  Heading=atan2(r_No_WP1_N_y, r_No_WP1_N_x); % Heading
%
% Tailcone Specific Constants:
% NOTE: Assume that for Yaw/Pitch control purposes, the ring and struts have 
% the properties of an Odyssey fin, with 4 times the surface area.  For 
% the surge equation, use the total (approximate, of course) area.
  xPcp= -3*in2m;            % Distance from pivot to ring CP in tailcone coords
  aR  = 2.50*in2m;          % Dist from c.l. to rudder actuator base
  rR1 = 1.00*in2m;          % See notes 99/4/20 p.13.  Values from Mark's 
  rR2 = 2.50*in2m;          % diagrams, p.16-17.                          
  xRR = 11.5*in2m;          % Center-hole to center-hole act. dist. in P.  zz
  maxDeltaR = pi/12;        % Max tailcone deflection; 15 degrees.        
% maxDeltaRDot = 3*pi/180;  % Max tailcone deflection rate; 5 deg/sec.    
% deltaRDead = .01*in2m/rR2;% Rudder deadband; 10 "mils"; actu. & pin.    
% deltaRDead = maxDeltaRDot*Ts;% Max rudder displacement in one sample interval
  maxActDisp = 1*in2m*m2c;  % One inch max displacement; counts
  maxActRate = .5*in2m*m2c; % Max actuator rate; counts/sec.
  maxActDispTs = maxActRate*Ts;% Max actuator disp. in one sample interval
  iSB  =0;                  % Current to power stepper controller board.
%
  StrutArea  = 0.0299;      % Odyssey fin area; m^2
  ControlArea = StrutArea*4;% Equivalent rudder/elevator area (m^2). 
  TotalArea = .47*pi*5*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 
%
% Sim Parameters:
  tStop = 2*60;              % Stop time, sec
  Npts  = 600;               % Number of points; should be tStop/Ts
%
% Open loop rudder commands:
  deltaR = 1*d2r;
  deltaR = 10*d2r;
  tdeltaR= 2;               % Rudder slew time
  tdur= 20;                 % Rudder slew time
%
% Control System:
  Kph = .7;                 % Unitless
  Krh = 2;                  % Sec; w/o integrator
  MaxRate = .2 * Krh;       % Rad/sec
%  Kwp = .0005;              % Rad/sec/meter; should be > 0; for Prop. Nav.
  Kwp = 0;                  % Rad/meter; should be > 0.
  Kwp = .05;                % Rad/meter; should be > 0.
  Kiwp= 0;                  % Rad/sec/meter; should be > 0. For integrator ver.
  Kiwp= .003;               % Rad/sec/meter; should be > 0. For integrator ver.
  WpMode = 0;               % 0=> Waypoint control; 1 => heading control; 
%
% Run the setup file here, if there is one.

% setHeading180
  setOpenLoop

%
% Computed Constants: Tailcone Actuators.
  disp0 = sqrt( (rR1 - xRR)^2 + (rR2 - aR)^2 ); %Act length at deltaR=0
  bAct  = sqrt(rR1^2 + rR2^2);   % radius of circle generated by pivot pt.
  d1Act = sqrt( aR^2 + xRR^2);   % Dist as shown on Act. notes, p. 21
  delta0= acos( (d1Act^2 + bAct^2 - disp0^2)/(2*d1Act*bAct) ); %See p.21
  ampsAct = .5/in2m*c2m*maxActRate+.9; %.5 is amps/(in/sec) from GT's diagram;
                                       %.9 is the y intercept.
%
% Open loop rudder profile
  tvals = [0, tdeltaR, tdeltaR+tdur, 2*tdeltaR+tdur, tStop];
  ovals = [0, deltaR,  deltaR,       0,              0 ];
% ovals = [0, 0,       0,            0,              0 ];
%
% Computed Constants: Vehicle Dynamics.
  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
  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
  Xudot = -.1*mass;         % Deriv of X wrto u-dot (Added Mass); kg; WAG
  Xudot = -1.583e-3*rho/2*len^3;% Deriv of X wrto u-dot; kg; DSRV
  Xrv   = .003358*rho/2*len^3;  % Deriv of X wrto rv; kg; DSRV
  Xrr   = -1.81e-4*rho/2*len^4; % Deriv of X wrto rr; N-sec^2; DSRV
  Xuu   = -5.03e-3*rho/2*len^2; % Deriv of X wrto uu; kg/m; DSRV
  Xuu   = -.5*rho*pi*(dia/2)^2*.15;  % Knut's computation from Brooks et. al.
  Xuu   = .5*rho*ControlArea*Cd0 - Tp/vel^2;  % See Notes, Eqn. 101.
  Xvv   = .0098*rho/2*len^2;    % Deriv of X wrto vv; kg/m; ODYSSEY
  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
%
% 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 = 1;                           % 1=> Plot linear as well; 0=> sim only
%
% Start screen log:
  eval(['!rm run', caseName,'.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:')
    tStop   = 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:')
    tStop  = 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(' -------------Open Loop Linear Analysis:-------------')
  disp(' ')
%
% Set up time vector:
  tsim = [0:Npts]'*tStop/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 = (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))
%
%
% -------------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;
              MSurge;
              vel;
              xR;
              xP;
              aR;
              rR1;
              rR2;
              xRR;
              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
              Xrv;
              Xrr;
              Xuu;
              Xvv;
              Yv;
              Yr;
              Yvabsr;
              Yvabsv;
              xG;
              Nv;
              Nr;
              Nrabsv;
              Nvabsv;
              Nrabsr;
              Decouple;
            ]';
%
% Run it!

% [tout] = sim('AltexSim1',tsim,simset('maxstep',tsim(2)-tsim(1)));
% [tout] = sim('AltexSim1',tsim,simset('maxstep',Ts));
  [tout] = sim('AltexSim2',tsim,simset('maxstep',Ts));
%
% -------------End Simulink--------------
%
%
% Indirect addressing for the columns of yout.
  lbls = { 'vIx'
           'rIx'
           'betaIx'
           'betaRIx'
           'AofAIx'
           'TrqNRIx'
           'TrqPIx'
           'xPosIx'
           'yPosIx'
           'yawIx'
           'FYRIx'
           'FYTIx'
           'FRAIx'
           'deltaRIx'
           'uIx'
           'rAmpsIx'
           'rAmpSecIx'
           'rActCntsIx'
           'rActDotIx'
         };
  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 rSS < 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.
%
  rSS = yout(Npts+1,rIx);
  Ap = (xR*Yvabsv - Nvabsv)/mass/vel^2;
  Bp = (xR*Yv - Nv)/mass/vel;
  Cp = ((xR*Yr - Nr)/mass - xR)*rSS/vel + Nrabsr/mass*(rSS/vel)^2;
  vSS = (-Bp + sqrt(Bp^2-4*Ap*Cp))/2/Ap;  %Positive root of a quadratic
  FYRSS = mass*vel*rSS - Yr*vel*rSS - Yv*vel*vSS - Yvabsv*vSS^2;
  turnDia = 2*vel/abs(rSS);
%
% Compute uSS, assuming that the Surge Eqn is uncoupled from the 2DOF.  That
% is, all of the U's (=uSS = vel) in the 2DOF are constant.
%
  vSS1 = yout(Npts+1,vIx);
  uSS = sqrt(-(mass*vSS1*rSS + Xvv*vSS1^2 + Xrv*vSS1*rSS + Xrr*rSS^2 + Tp)/...
              (Xuu - .5*rho*ControlArea*Cd0/pi/Ar/Ef));
  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.', rSS/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*rSS)/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)));
  disp(sprintf(' %-15f %-15f: uSS [Sim Comp], m/sec.',yout(Npts+1,uIx),uSS));
%
  diary off    
  clear functions mex

  plotAltex

  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


