close all
clear all

hard_copy = 0;

data = load('2-5May11.CSV');
data_b= load('2b-5May11.CSV');

t = [data(:,1); data(end,1)+.02+data_b(:,1)];
x_raw = [data(:,2); data_b(:,2)];
F_raw = [data(:,3); data_b(:,3)];
MP1_raw = [data(:,4); data_b(:,4)];
MP2_raw = [data(:,5); data_b(:,5)];
V_raw = [data(:,6); data_b(:,6)];
LCP_raw = [data(:,7); data_b(:,7)];
I_raw = [data(:,8); data_b(:,8)];
N_raw = [data(:,9); data_b(:,9)];

date = [data(:,10); data_b(:,10)];
time = [data(:,11); data_b(:,11)];

clear data data_b;

x_raw = (x_raw-10.42)*59.4/60.92;


windowSize = 3;
x = filtfilt(ones(1,windowSize)/windowSize,1,x_raw);
xdot = filtfilt(ones(1,windowSize)/windowSize,1,[0; diff(x)./diff(t)]);


figure('Position',[ 21         273        1536         820],'PaperOrientation','landscape','PaperPosition',[0 0 11 8.5]);
plot(xdot);
ylabel('ram speed (ips)');
xlabel('measurement step');
foo = axis;
foo(3) = -10;
foo(4) = 20;
axis(foo);
if (hard_copy == 1)
  print('-dpdf','RawSpeedHistory');
  print('-dpng','RawSpeedHistory');
end



F = F_raw;
MP1 = MP1_raw;
MP2 = MP2_raw;
LCP = LCP_raw;
V = V_raw;
I = I_raw;
%MP1 = filtfilt(ones(1,windowSize)/windowSize,1,MP1_raw);
%MP2 = filtfilt(ones(1,windowSize)/windowSize,1,MP2_raw);
%LCP = filtfilt(ones(1,windowSize)/windowSize,1,LCP_raw);
%I = filtfilt(ones(1,windowSize)/windowSize,1,I_raw);
%V = filtfilt(ones(1,windowSize)/windowSize,1,V_raw);
%N = filtfilt(ones(1,windowSize)/windowSize,1,N_raw);


Ap = .25*pi*(4-1.375^2);

x = x*.0254;  %meters
xdot = xdot*.0254;  %m/s
Ap = Ap*(.0254^2);
F = F*4.45;

MP1 = MP1*6894;
MP2 = MP2*6894;

resistances = [inf 50 25 10 5 3.5 0]; %ohms
speeds = [1:18];  %ips
n_loads = length(resistances);
n_speeds = length(speeds);
idx = cell(n_speeds,n_loads);

%speed = 1ips
idx{1,2} = 68800:71680;
idx{1,3} = 73010:75830;
idx{1,4} = 76710:79550;
idx{1,5} = 80310:83170;
idx{1,6} = 84420:87220;
idx{1,1} = 122020:124890;

%speed = 2ips
idx{2,2} = 537:1726;
idx{2,3} = 2884:4139;
idx{2,4} = 5167:6406;
idx{2,5} = 7526:8700;
idx{2,6} = 10050:11320;
idx{2,1} = 125780:127090;


%speed = 3ips
idx{3,2} = 89320:90240;
idx{3,3} = 91810:92690;
idx{3,4} = 94280:95180;
idx{3,5} = 96500:97390;
idx{3,6} = 99430:100350;
idx{3,1} = 127980:128890;

%speed = 4ips
idx{4,2} = 12850:13420;
idx{4,3} = 14530:15050;
idx{4,4} = 16590:17190;
idx{4,5} = 18140:18750;
idx{4,6} = 20000:20560;
idx{4,1} = 129630:130280;

%speed = 5ips
idx{5,2} = 102309:102810;
idx{5,3} = 105064:105595;
idx{5,4} = 106982:107510;
idx{5,5} = 109009:109539;
idx{5,6} = 111130:111650;
idx{5,1} = 131030:131540;

%speed = 6ips
idx{6,2} = 21790:22170;
idx{6,3} = 23420:23770;
idx{6,4} = 24670:25040;
idx{6,5} = 26090:26440;
idx{6,6} = 27750:28110;
idx{6,1} = 132280:132680;

%speed = 7ips
idx{7,2} = 113170:113540;
idx{7,3} = 114730:115090;
idx{7,4} = 116270:116640;
idx{7,5} = 118030:118390;
idx{7,6} = 120030:120400;
idx{7,1} = 133670:133990;

%speed = 8ips
idx{8,2} = 29840:30140;
idx{8,3} = 31040:31340;
idx{8,4} = 32220:32530;
idx{8,5} = 33630:33930;
idx{8,6} = 34900:35190;
idx{8,1} = 134960:135260;

%speed = 10ips
idx{10,2} = 36290:36510;
idx{10,3} = 37760:37990;
idx{10,4} = 38880:39120;
idx{10,5} = 40140:40380;
idx{10,6} = 44190:44410;
idx{10,1} = 136800:137020;

%speed = 12ips
idx{12,2} = 45900:46080;
idx{12,3} = 47190:47360;
idx{12,4} = 48410:48590;
idx{12,5} = 49630:49820;
%idx{12,6} = :;
idx{12,1} = 138090:138260;

%speed = 14ips
idx{14,2} = 51030:51190;
idx{14,3} = 52340:52500;
idx{14,4} = 53860:54020;
idx{14,5} = 55250:55390;
%idx{14,6} = :;
idx{14,1} = 139210:139330;

%speed = 16ips
idx{16,2} = 56800:56920;
idx{16,3} = 58220:58340;
idx{16,4} = 59330:59450;
idx{16,5} = 61100:61230;
%idx{16,6} = :;
idx{16,1} = 140410:140510;

%speed = 18ips
idx{18,2} = [63320:63430 64590:64690];  %This run was repeated
idx{18,3} = 66200:66300;
idx{18,4} = 67300:67400;
%idx{18,6} = :;
idx{18,1} = 141710:141770;

%Electrical Frequency Data
elec_f = NaN(18,6);

%speed = 1 ips
elec_f(1,2) = 22; %Hz
elec_f(1,3) = 22;
elec_f(1,4) = 21;
elec_f(1,5) = 22;
elec_f(1,6) = 22;
elec_f(1,1) = 22;

%speed = 2 ips
elec_f(2,2) = 44;  %Hz
elec_f(2,3) = 44;
elec_f(2,4) = 44;
elec_f(2,5) = 44;
elec_f(2,6) = 43;
elec_f(2,1) = 45;

%speed = 3 ips
elec_f(3,2) = 67; %Hz
elec_f(3,3) = 66;
elec_f(3,4) = 66;
elec_f(3,5) = 66;
elec_f(3,6) = 66;
elec_f(3,1) = 68;

%speed = 4 ips
elec_f(4,2) = 88; %Hz
elec_f(4,3) = 87;
elec_f(4,4) = 88;
elec_f(4,5) = 87;
elec_f(4,6) = 87;
elec_f(4,1) = 88;

%speed = 5 ips
elec_f(5,2) = 110; %Hz
elec_f(8,3) = 110;
elec_f(5,4) = 108;
elec_f(5,5) = 109;
elec_f(5,6) = 108;
elec_f(5,1) = 110;

%speed = 6 ips
elec_f(6,2) = 134; %Hz
elec_f(6,3) = 133;
elec_f(6,4) = 131;
elec_f(6,5) = 130;
elec_f(6,6) = 129;
elec_f(6,1) = 132;

%speed = 7 ips
elec_f(7,2) = 155; %Hz
elec_f(7,3) = 154;
elec_f(7,4) = 153;
elec_f(7,5) = 153;
elec_f(7,6) = 152;
elec_f(7,1) = 159;

%speed = 8 ips
elec_f(8,2) = 176; %Hz
elec_f(8,3) = 177;
elec_f(8,4) = 175;
elec_f(8,5) = 175;
elec_f(8,6) = 175;
elec_f(8,1) = 178;

%speed = 10 ips
elec_f(10,2) = 220; %Hz
elec_f(10,3) = 221;
elec_f(10,4) = 219;
elec_f(10,5) = 217;
elec_f(10,6) = 219;
elec_f(10,1) = 218;

%speed = 12 ips
elec_f(12,2) = 262; %Hz
elec_f(12,3) = 262;
elec_f(12,4) = 262;
elec_f(12,5) = 261;
%elec_f(12,6) = ;
elec_f(12,1) = 261;

%speed = 14 ips
elec_f(14,2) = 306; %Hz
elec_f(14,3) = 307;
elec_f(14,4) = 306;
elec_f(14,5) = 306;
%elec_f(14,6) = ;
elec_f(14,1) = 311;

%speed = 16 ips
elec_f(16,2) = 351; %Hz
elec_f(16,3) = 350;
elec_f(16,4) = 349;
elec_f(16,5) = 351;
%elec_f(16,6) = ;
elec_f(16,1) = 351;

%speed = 18 ips
elec_f(18,2) = 395; %Hz
elec_f(18,3) = 390;
elec_f(18,4) = 391;
%elec_f(18,5) = ;
%elec_f(18,6) = ;
elec_f(18,1) = 398;

N = elec_f*15;  %RPM

for i = 1:18
for j = 1:6
  force(i,j) = mean(F(idx{i,j}));
  speed(i,j) = mean(xdot(idx{i,j}));
  damping(i,j) = mean(F(idx{i,j})./xdot(idx{i,j}));
  sr_damping(i,j) = mean(F(idx{i,j})./(xdot(idx{i,j}).^.5));
  mot_pres(i,j) = mean(MP1(idx{i,j})-MP2(idx{i,j}));
  mech_pow(i,j) = mean(F(idx{i,j}).*xdot(idx{i,j}));
  elec_pow(i,j) = mean(V(idx{i,j}).*I(idx{i,j}));
  hyd_pow(i,j) = mean(Ap*xdot(idx{i,j}).*(MP1(idx{i,j})-MP2(idx{i,j})))
  Resistance(i,j) = mean(V(idx{i,j})./I(idx{i,j}));
  N_dataq(i,j) = mean(N_raw(idx{i,j}));
  
  eff(i,j) = elec_pow(i,j)./mech_pow(i,j);
  hyd_eff(i,j) = hyd_pow(i,j)./mech_pow(i,j);
  motgen_eff(i,j) = elec_pow(i,j)./hyd_pow(i,j);
end
end

figure('Position',[520   211   767   887],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
plot(speeds,speed*39.4,'kx')
xlabel('Specified Speed (ips)');
ylabel('Measured Speed (ips)');
grid on
title('Test Machine Performance')

subplot(2,1,2)
plot(resistances(2:6),Resistance(:,2:6),'kx')
xlabel('Selected Resistor Value');
ylabel('V/I');
grid on
if (hard_copy == 1)
  print('-dpdf','TestMachinePerformance');
  print('-dpng','TestMachinePerformance');
end



figure('Position',[520   211   767   887],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
plot(N([1:8 10:2:18]),N([1:8 10:2:18])*.99,N([1:8 10:2:18],:),N_dataq([1:8 10:2:18],:));
xlabel('RPM from oscilloscope');
ylabel('RPM from Dataq (averaged over run)');
title('RPM Measurement Accuracy')
legend('99.0% derate','inf','50 Ohms','25 Ohms','10 Ohms','5 Ohms','3.5 Ohms','Location','NorthWest');
if (hard_copy == 1)
  print('-dpdf','RPMMeasurementAccuracy');
  print('-dpng','RPMMeasurementAccuracy');
end


no_leak_rpms = Ap*speed*60/(4.9/1000000);
eff_v = .99*N./no_leak_rpms;

figure('Position',[520   211   767   887],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
plot(speed([1:8 10:2:18])/.0254,no_leak_rpms([1:8 10:2:18]),'k-',speed([1:8 10:2:18],:)/.0254,.99*N([1:8 10:2:18],:))
xlabel('Measured Ram Speed (ips)');
ylabel('Measured Shaft Speed (rpm)');
legend('No Leak','inf','50 Ohms','25 Ohms','10 Ohms','5 Ohms','3.5 Ohms','Location','NorthWest');
subplot(2,1,2)
plot(speed([1:8 10:2:18],:)/.0254,eff_v([1:8 10:2:18]));
xlabel('Measured Ram Speed (ips)');
ylabel('Measured Volumetric Efficiency');



figure('Position',[520   211   767   887],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
plot(speed([1:8 10:2:18],:),force([1:8 10:2:18],:),speed([1:8 10:2:18]),20000*(.0254*speeds([1:8 10:2:18])).^.5,'--',speed([1:8 10:2:18]),14000*(speed([1:8 10:2:18])).^.5,'--',speed([1:8 10:2:18]),900+7000*(speed([1:8 10:2:18])),'--');
xlabel('speed (m/s)');
ylabel('Force (N)');
title('System Force and Efficiency');
hleg = legend('inf','50 Ohms','25 Ohms','10 Ohms','5 Ohms','3.5 Ohms','20000 N/sqrt(m/s)','14000 N/sqrt(m/s)','900N+7000 N/(m/s)','Location','NorthWest');
set(hleg,'FontSize',8)
subplot(2,1,2)
plot(speed([1:8 10:2:18],:),eff([1:8 10:2:18],:))
ylabel('System Efficiency')
if (hard_copy == 1)
  print('-dpdf','ForceAndEfficiency');
  print('-dpng','ForceAndEfficiency');
end


figure('Position',[520   211   767   887],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
plot(speed([1:8 10:2:18],:),1.66*mot_pres([1:8 10:2:18],:)/6894-force([1:8 10:2:18],:)/4.45);

legend('inf Ohms','50 Ohms','25 Ohms','10 Ohms','5 Ohms','3.5 Ohms');
xlabel('speed (m/s)');
ylabel('Seal and Fluid Drag (lbs)(A*deltaP-F)');
subplot(2,1,2)
plot(speed([1:8 10:2:18],:),hyd_eff([1:8 10:2:18],:))
ylabel('Plumbing Efficiency')
xlabel('speed (m/s)');

if (hard_copy == 1)
  print('-dpdf','PlumbingEfficiency');
  
end

figure('Position',[520   211   767   887],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
plot(speed([1:8 10:2:18],:),motgen_eff([1:8 10:2:18],:))
legend('inf Ohms','50 Ohms','25 Ohms','10 Ohms','5 Ohms','3.5 Ohms');
ylabel('Motor/Generator Efficiency')
xlabel('m/s');
if (hard_copy == 1)
  print('-dpdf','MotorGenEfficiency');
end