%
% PURPOSE: Plot tank test data in Matlab.  Compare Semi and Moog
% AUTHOR:  Rob McEwen
% DATE:  98/12/15
% MODS: 
% 98/11/19:
% 1) Divided V2lbs by 2 to account for Gary's voltage divider.
% 2) The conditioning amp zero is now at 2.5 volts, or 2048 A/D counts.
% Full scale is 0 to 5 volts. Previously, the cond. amp output range 
% was 106.38 lbs/Volt, centered on about -300 mV.
%
% 
  nameStr = input(' Data File Name: ','s')
%
% Set this to true to configure for the SemiPower-to-SemiPower comparison.
  S2S = 1;     
%
% Load the data
  if(nameStr) 
    if (~S2S)
      eval([nameStr,'_VEL10'])        %Moog w/ VEL control, 100 Hz.
    else
      eval([nameStr,'_SemiVEL50_2'])  %SemiPower w/ VEL cntrl @ 20Hz, long durn
    end
    nameStrLast = nameStr;
    doit = 1;
  else
    nameStr = nameStrLast;
    doit = 0;
  end
%
% Add lables below so that they match the order of the columns in asciiData.
%
% cnts2V = 5/2^11;           %Datel Conversion from AD counts to volts.
  cnts2V = 5/2^12;           %B&B Conversion from AD counts to volts.
  V2lbs  = 1/.00940*2;       %Conversion from Volts to lbs, from cal. data.
  lbs2N  = 4.448;            %Newtons per pound.
%
  c2N    = cnts2V*V2lbs*lbs2N;
  v2N    = V2lbs*lbs2N;

  lbls = { 'time'
           'ADstatus'
           'ADvalue'
           'vel'
           'temp'
           'thrust'
           'cmdThr'
           'torque'
         };
  sz = size(lbls);
  i = 1;
  while( i <= sz(1) ),
    eval([ char(lbls(i))  ' = '  num2str(i) ';']);
    i = i + 1;
  end
%
% Load the SemiPower data:
  if(doit) 
    dataMoog = data;
    if (~S2S)
      eval([nameStr,'_VEL50'])      %SemiPower w/ VEL control, short duration
    else
      eval([nameStr,'_SpeedReg'])   %SemiPower w/ OnBoard speed regulator
    end
    dataSemi = data;
    clear data;
  end
%
% First, compute the initial (first 10 seconds) strain gauge bias reading:
  [mm,idxS] = min(abs(dataSemi(:,time) - 10 ));
  biasS = ones(1,idxS)*dataSemi(1:idxS,ADvalue)/idxS*c2N;
  sgaugeS = dataSemi(:,ADvalue)*c2N - biasS;
%
  [mm,idxM] = min(abs(dataMoog(:,time) - 10 ));
  biasM = ones(1,idxM)*dataMoog(1:idxM,ADvalue)/idxM*c2N;
  sgaugeM = dataMoog(:,ADvalue)*c2N - biasM;
%
  disp([' The (Semi Moog) bias is ', num2str([biasS biasM],4), ' N, rms.'])
%
% Now, compute the RMS noise in the first 10 seconds.
  NoiseS = std(sgaugeS(1:idxS));
  NoiseM = std(sgaugeM(1:idxM));
  disp(' ')
  disp([' The (S,M) noise is ', num2str([NoiseS NoiseM],4), ' Newtons, rms.'])
%
% Data strings in the plot:  
  textStr = strvcat(myDate,...
            ['The bias is ',num2str(biasS,4),' ',num2str(biasM,4),' N'],...
            ['The 0-10 second noise is ',num2str(NoiseS,3),' '],...
            [num2str(NoiseM,3),' N, rms'],...
            ['File Name: ',nameStr]);
%
% Reconcile different sampling rates between the Semi and Moog events.
  szS = size(dataSemi);
  szM = size(dataMoog);
  if szS(1) == szM(1)
     dataSemiD = dataSemi;
     sgaugeSD  = sgaugeS;
  else
     %Check to see if sampling rate was 50 Hz, vs 100 Hz.
     if szS(1) == 2*szM(1)
        dataSemiD = dataSemi(1:2:szS(1),:);
        sgaugeSD  = sgaugeS(1:2:szS(1));
     else
        disp(' Array Dimension Error.')
     end
  end
%
% Set up plot label strings
  xlabelStr = ['Seconds.  10 Sec Noise (S M) = ', ...
               num2str([NoiseS NoiseM],3),' N, rms.'];
  titleStr = [myDate,' ',nameStr,' (S M) Bias: ', ...
  num2str([biasS biasM],4),' Cnts'];

  Np = 3;                          %Number of plots
  figure(1)
  clf
  sz = size(dataSemiD);
  N = sz(1);
  plot(dataSemiD(1:N,time), dataSemiD(1:N,vel),...
       dataSemiD(1:N,time), dataMoog(1:N,vel),'-')
  ylabel('Shaft Speed, RPM')
  xlabel(xlabelStr)
  title(titleStr)
  grid on
% val = axis;
% text(aval(2)*.025,aval(4)*.90,textStr)
  if (~S2S) 
    legend('Semi','Moog',0)
  else
    legend('Speed Reg','Vel 20Hz',0)
  end

  figure(2)
  clf
  sz = size(dataSemiD);
  N = sz(1);
  plot(dataSemiD(1:N,time), sgaugeSD(1:N),...
       dataSemiD(1:N,time), sgaugeM(1:N),'-',...
       dataSemiD(1:N,time), dataSemiD(1:N,thrust),'-')
  ylabel('Thrust, Newtons')
  xlabel(xlabelStr)
  title(titleStr)
  grid on
% val = axis;
% text(aval(2)*.025,aval(4)*.90,textStr)
  if (~S2S) 
    legend('Semi','Moog','Commanded',0)
  else
    legend('Speed Reg','Vel 20Hz','Commanded',0)
  end

  figure(3)
  clf
  sz = size(dataSemiD);
  N = sz(1);
  plot(dataSemiD(1:N,time),[dataSemiD(1:N,torque)/100])
  ylabel('SemiPower Torque, Newton-meters')
  xlabel(['Seconds.  10 Sec Noise = ', num2str(NoiseS), ' Newtons, rms.'])
  title([myDate,' ',nameStr,' Bias: ', num2str(biasS), ' Counts'])
  grid on
%
% Plot a smaller portion:
%
%  eachChar = (nameStr(1:5) ~= 'SlowT') & (nameStr(1:5) ~= 'NegSt' );
   eachChar = (nameStr(1:5) ~= 'NegSt' );
%
%  eachChar is 1x5 and has a one for each char that is not equal.
%
   if ( (eachChar * ones(5,1)) > 0 )

    Np = 6;
    if nameStr(1:5) == 'IncTr', N1=751; N2=1001; end
    if nameStr(1:5) == ('IncSq' | 'SlowT' ), N1=6001;N2=7001; end
    if nameStr(1:5) == 'Steps', N1=24901;N2=25201; end
    
    if (S2S)
      if nameStr(1:5) == 'IncTr', N1=4500; N2=5001; end
      if (nameStr(1:5) == 'IncSq') |(nameStr(1:5) ==  'SlowT' ), 
        N1=6001;N2=7001; end
      if nameStr(1:5) == 'Steps', N1=24901;N2=25201; end
    end   

    
    figure(4)
    clf
    plot(dataSemiD(N1:N2,time), dataSemiD(N1:N2,vel),...
    dataSemiD(N1:N2,time), dataMoog(N1:N2,vel),'-')
    ylabel('Shaft Speed, RPM')
    xlabel(xlabelStr)
    title(titleStr)
    grid on
    % val = axis;
    % text(aval(2)*.025,aval(4)*.90,textStr)
    if (~S2S) 
      legend('Semi','Moog',0)
    else
      legend('Speed Reg','Vel 20Hz',0)
    end

    figure(5)
    clf
    plot(dataSemiD(N1:N2,time), sgaugeSD(N1:N2),...
    dataSemiD(N1:N2,time), sgaugeM(N1:N2),'-',...
    dataSemiD(N1:N2,time), dataSemiD(N1:N2,thrust),'-')
    ylabel('Thrust, Newtons')
    xlabel(xlabelStr)
    title(titleStr)
    grid on
    % val = axis;
    % text(aval(2)*.025,aval(4)*.90,textStr)
    if (~S2S) 
      legend('Semi','Moog','Commanded',0)
    else
      legend('Speed Reg','Vel 20Hz','Commanded',0)
    end

    figure(6)
    clf
    plot(dataSemiD(N1:N2,time),[dataSemiD(N1:N2,torque)/100])
    ylabel('SemiPower Torque, Newton-meters')
    xlabel(['Seconds.  10 Sec Noise = ', num2str(NoiseS), ' Newtons, rms.'])
    title([myDate,' ',nameStr,' Bias: ', num2str(biasS), ' Counts'])
    grid on
   end

  f7exist = 0;
  szStr = size(nameStr);
  if szStr(2) >= 6
     nameStrR = nameStr(1:6);
     if nameStrR == 'SlowTr' | nameStrR == 'IncTri'
       Np = 7;
       figure(7)
       clf
       sz = size(dataSemiD);
       N = sz(1);
       plot(dataSemiD(1:N,vel).*abs(dataSemiD(1:N,vel)), sgaugeSD,...
            dataMoog(1:N,vel).*abs(dataMoog(1:N,vel)),  sgaugeM,'g:')
       ylabel('Torque, Newton-meters')
       xlabel([' RPM^2 '])
       grid on
       title(titleStr)
       if (~S2S) 
	 legend('Semi','Moog',0)
       else
	 legend('Speed Reg','Vel 20Hz',0)
       end
       f7exist = 1;
     end
  end

  doit = input(' Hardcopy? y => .eps, Y => printer ','s');
  if(doit == 'y' | doit == 'Y')
    for i = 1:Np
      eval(['print -f',num2str(i),'  -depsc ',nameStr,'_',num2str(i)])
    end
    if f7exist, eval(['print -f7 -depsc ',nameStr,'_7']), end
    if( doit == 'Y' )
      for i = 1:Np
	eval(['!lpr ',nameStr,'_',num2str(i),'.eps'])
      end
      if f7exist, eval(['!lpr ',nameStr,'_7.eps']), end
    end
  end

return

  figure(1)
  clf

  subplot(3,1,1), 
  plot(data(:,time),data(:,thrust))
  ylabel('Newtons')
  title(['Cmded Thrust.',' File: ',nameStr,' Bias: ', num2str(bias),' ',myDate])

  subplot(3,1,2), 
  plot(data(:,time),data(:,vel))
  ylabel('RPM')
  title('Propeller Speed')

  subplot(3,1,3), 
  plot(data(:,time),data(:,temp)),
  xlabel('Time, seconds')
  ylabel('Deg C')
  title('Motor Temperature')

%
% Filter that baby:
  [BB,AA] = butter(4,1/3);  %15/3 = 5 Hz cutoff freq; 4th order
  sgf = filter(BB,AA,sgauge);

  figure(3)
  clf
  sz = size(data);
  N = sz(1);
  plot(data(1:N,time),[data(1:N,thrust) sgf(1:N)])
  ylabel('Thrust, Newtons; 5Hz 4pole Butterworth')
  xlabel([' Seconds'])
  grid on
  legend('Commanded','Measured')
  title([myDate,' File: ',nameStr,' Bias: ', num2str(bias)])
% text(aval(2)*.025,aval(4)*.90,strvcat(textStr,'5 Hz 4th Order Butterworth'))

  figure(3)
  clf

  subplot(2,1,1), 
  plot(data(:,time),sgauge );
  ylabel('Newtons')
  title('Strain Gauge')

  subplot(2,1,2), 
  plot(data(:,time),data(:,ADstatus))
  ylabel('Counts')
  title('Datel A/D Board Status')
  xlabel('Time, seconds')
