%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PURPOSE: Plot i versus U for the two thruster curves from the Ody sim.
%          Test reverse
% AUTHOR:  Rob McEwen
% DATE:    00/1/4
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   /*------------------------------------------------------------------*
%     The old model uses a drag power that is quadratic in desired speed,
%     which contradicts hydrodynamic theory. This model uses a cubic power
%     law that combines the fudge into effective C_d.
%
%                      I_des = 1/2 rho U^3 Cd A / V_batt
%
%     V_batt ~ 47.5 under load, A ~ 0.25 m^2, rho = 1025, combine all
%     constants yields 2.6978
%    *------------------------------------------------------------------*/
%

  clear U Xuabu LHS

  V_batt  = 47.5;      %Volts
  Cdf     = .21;       %Forward drag coeff; unitless
  rho     = 1025;      %Density of Seawater; kg/m^3
  Af      = .25;       %Frontal area; m^2
  Kt      = .456;      %torque constant (in Nm/A) - inc gearbox */
  Eta	  = 0.80;      %propellor efficiency = KT/KE */
  Pit     = 0.124;     %prop pitch: advance/rad rotation (m/rad)*/
  Area    = 0.28;      %water jet (duct) projected area (in m2) */

%
% Computed constants:
%
  Xuabu   = -.5*rho*Cdf*Af;
  Xuabu   = -9.1;
  U       = [0:500]'/500*(1.75-1.25) + 1.25; % Speed vector - indep axis; m/s
  U       = [0:500]'/500*3 -1.5; % Speed vector - indep axis; [-1.5, 1.5] m/s

%
% thruster_current  = 2.6978 * sensor(EFFECTIVE_DRAG_COEFF) * pow(des_speed,3);
%
  thruster_current1 = .5*rho*Af*Cdf* U.^3/V_batt;
%
  LHS = -Xuabu*(U.*abs(U)).*( 2.455 - 2.175*U +  0.75*(U.*abs(U)));
  LHS = -Xuabu*(U.^2).*( 2.455 - 2.175*U +  0.75*(U.^2));

  des_omega = ( Area*rho*Eta*Pit*U + ...
                sqrt( (Area*rho*Eta*Pit*U).*abs(Area*rho*Eta*Pit*U) + ...
                       4.*LHS*Area*rho*Eta*Eta*Pit*Pit ) ) / ...
                      (2*Area*rho*(Eta*Pit)^2);

  thruster_current2 = (.5*rho*Area*(Eta*Pit)^3.*des_omega.*abs(des_omega) - ...
                            .5*Area*rho*Eta*Pit*U.*abs(U))/Kt.*sign(U);
  %/* make sign same as desired speed */
  %thruster_current2 *= sgn(U);
  

  caseName = 'QNX';
  str = [', Case:',caseName,' RSM ',myDate];
  PaperPosn = [.5 .5 3.75 2.5];
  FigurePosn = [5 30 815 760];
  i = 0;

  i=i+1;
  figure(i)
  clf
  set(gcf,'position',FigurePosn)
  plot(U, [thruster_current1 thruster_current2 des_omega]);
  grid on
  legend('i Now','i Theory','omega, Theory')
  ylabel('Amperes & Rad/Sec')
  xlabel('Speed, Meters/Second')
  title(['Thruster Current Versus Speed',str]);

  nameStr = 'thrCompPlt';
  doit = input(' Hardcopy? y => Create .eps, Y => Also send to Tahoe ','s');
  if(doit == 'y' | doit == 'Y')
    eval(['print -f1 -depsc ',nameStr,'_1'])
    if( doit == 'Y' )
      eval(['!lp -opostscript -onb -dtahoe ',nameStr,'_1.eps'])
    end
  end


