close all
clear all

rho = 700;  %kg/m^3
nu = 2.5*15e-6;   %ISO-15 oil at 20 degrees C (~2.5 times the viscosity at 40 degrees C).

pipe_size = [1 2];  %Nominal pipe sizes
pipe_id = [.584 .620];
K_elbow = [0.69 0.57];
K_Tee_FlowThrough = [0.46 0.38];
K_Tee_FlowBranched = [1.38 1.14];
Cv_CheckValve = [12 60];
K_ent = 0.5;
K_exit = 1.0;


pipe_size = 1;

D = 2.5*.0254;
Dr = 2*.0254;
d = pipe_id(pipe_size)*0.0254;
e = 1.5E-6;  %m
L = 2;  %m


A = 0.25*pi*(D^2-Dr^2);
Ap = 0.25*pi*d^2;

V = 0:1/100:1;
F = 0:10:15000; F = F';


b = 14000;
figure
plot(V,(b*sqrV/A)/6894);  %psi
figure
plot(V,V*A*(39.4^3));


PressIdeal = F/A;
Q = V*A;

PowIdeal = PressIdeal*Q;

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
surf(V,F,PowIdeal,'LineStyle','none')
xlabel('Ram Speed (m/s)');
ylabel('Force (N)');
zlabel('Watts');
title('Input Power as a function of ram speed and force');
colorbar


%Compute headloss for non-rectified case  - Entire Force-Velocity Space
Vp = Q/Ap;

Re = V*D/nu;
Rep = Vp*d/nu;


f1 = friction_factor(Re,e/D);
f2 = friction_factor(Rep,e/D);


hl(1,:) = L*f1.*V.^2/(2*D);
hl(2,:) = (K_ent*Vp.^2)/2;
hl(3,:) = (K_elbow(pipe_size)*Vp.^2)/2;
hl(4,:) = L*f2.*Vp.^2/(2*d);
hl(5,:) = (K_elbow(pipe_size)*Vp.^2)/2;
hl(6,:) = (K_exit*Vp.^2)/2;


deltaP = rho*(sum(hl,1));
F_parasitic = deltaP*A;

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
subplot(2,1,1)
plot(V,rho*hl/6894.757,V,deltaP/6894.757)
legend('dP-Ram','dP-Entrance','dP-Elbow','dP-Pipe','dP-Elbow','dP-Exit','Total Pressure Drop','Location','NorthWest');
xlabel('Ram Speed (m/s)');
ylabel('Plumbing Pressure Drop (psi)');
title(['Losses for Non-Rectified Case - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);


subplot(2,1,2)
plot(V,hl*rho*A/4.45,V,F_parasitic/4.45);
legend('dF-Ram','dF-Entrance','dF-Elbow','dF-Pipe','dF-Elbow','dF-Exit','Total Force Loss','Location','NorthWest');
xlabel('Ram Speed (m/s)');
ylabel('Force (lbs)');



for i = 1:length(F)
  for j = 1:length(V);  
    PressAtMotor(i,j) = PressIdeal(i)-deltaP(j);
    PowAtMotor(i,j) = PressAtMotor(i,j)*Q(j);
  end
end

eff = PowAtMotor./PowIdeal;

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
surf(V,F,eff,'LineStyle','none')
zlabel('Efficiency');
xlabel('Ram Speed (m/s)');
ylabel('Input force (N)');
title(['Efficiences for Non-Rectified Case at All Damping Relationships - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);


%Compute losses and efficiency for specific relationships of Force and
%Velocity (slices through the 3-d plots.

%square root cases
clear eff PressAtMotor PowAtMotor PressIdeal PowIdeal F;
b = 8000:2000:16000; b = b'; %N/(m/s)^2
F = b*V.^.5;

PressIdeal = F/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);
eff = PowAtMotor./PowIdeal;
figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');

ah1 = subplot(3,1,1)
plot(V,eff);
legend('b = 8000 N/(m/s)^.5','b = 10000 N/(m/s)^.5','b = 12000 N/(m/s)^.5','b = 14000 N/(m/s)^.5','b = 16000 N/(m/s)^.5');
ylabel('Efficiency');

title(['Efficiences for Non-Rectified Case at Various Damping Relationships - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);

%Linear cases
clear eff PressAtMotor PowAtMotor;


PressIdeal = b*V/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);

eff = PowAtMotor./PowIdeal;
ah2 = subplot(3,1,2)
plot(V,eff);
legend('b = 8000 N/(m/s)','b = 10000 N/(m/s)','b = 12000 N/(m/s)','b = 14000 N/(m/s)','b = 16000 N/(m/s)');
ylabel('Efficiency');

%Quadratic cases
clear eff PressAtMotor PowAtMotor PressIdeal PowIdeal F;

F = b*V.^2;

PressIdeal = F/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);
eff = PowAtMotor./PowIdeal;
ah3 = subplot(3,1,3)
plot(V,eff);
legend('b = 8000 N/(m/s)^2','b = 10000 N/(m/s)^2','b = 12000 N/(m/s)^2','b = 14000 N/(m/s)^2','b = 16000 N/(m/s)^2');
ylabel('Efficiency');

ylabel('Efficiency');
xlabel('Ram Speed (m/s)');

linkaxes([ah1 ah2 ah3],'x');



%Compute headloss for rectified case - Entire Force-Velocity Space
F = 0:10:15000; F = F';

PressIdeal = F/A;
Q = V*A;

PowIdeal = PressIdeal*Q;

hl(1,:) = L*f1.*V.^2/(2*D);
hl(2,:) = (K_ent*Vp.^2)/2;
hl(3,:) = (K_Tee_FlowThrough(pipe_size)*Vp.^2)/2;
hl(4,:) = (K_elbow(pipe_size)*Vp.^2)/2;
hl(5,:) = 3*(K_Tee_FlowBranched(pipe_size)*Vp.^2)/2;
hl(6,:) = L*f2.*Vp.^2/(2*d);
hl(7,:) = (K_exit*Vp.^2)/2;
hl(8,:) = 2*(rho/1000)*(Q/Cv_CheckValve(pipe_size)).^2;


deltaP = rho*(sum(hl,1));
F_parasitic = deltaP*A;

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
subplot(2,1,1)
plot(V,rho*hl/6894.757,V,deltaP/6894.757)
legend('dP-Ram','dP-Entrance','dP-T_Straight (1)','dP-Elbow','dP-T-Branched (3)','dP-Pipe','dP-Exit','dP-CheckValves (2)','Total Pressure Drop','Location','NorthWest');
xlabel('Ram Speed (m/s)');
ylabel('Plumbing Pressure Drop (psi)');
title(['Losses for Rectified Case - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);


subplot(2,1,2)
plot(V,hl*rho*A/4.45,V,F_parasitic/4.45);
legend('dF-Ram','dF-Entrance','dF-T_Straight (1)','dF-Elbow','dF-T-Branched (3)','dF-Pipe','dF-Exit','dF-CheckValves (2)','Total Force Loss','Location','NorthWest');
xlabel('Ram Speed (m/s)');
ylabel('Force (lbs)');

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
plot(V,V*39.4/12,V,Vp*39.4/12);
legend('Ram flow speed','Pipe flow speed');
xlabel('Ram Speed (m/s)');
ylabel('Flow speed in ft/s')
title(['Average Flow Speeds - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);


for i = 1:length(F)
  for j = 1:length(V);  
    PressAtMotor(i,j) = PressIdeal(i)-deltaP(j);
    PowAtMotor(i,j) = PressAtMotor(i,j)*Q(j);
  end
end

eff = PowAtMotor./PowIdeal;
figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
surf(V,F,eff,'LineStyle','none')
zlabel('Efficiency');
xlabel('Ram Speed (m/s)');
ylabel('Input force (N)');
title(['Efficiences for Rectified Case at All Damping Relationships - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);


%Compute losses and efficiency for specific relationships of Force and
%Velocity (slices through the 3-d plots.

%square root cases
clear eff PressAtMotor PowAtMotor PressIdeal PowIdeal F;
b = 8000:2000:16000; b = b'; %N/(m/s)^2

F = b*V.^.5;

PressIdeal = F/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);
eff = PowAtMotor./PowIdeal;
figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
ah1 = subplot(3,1,1)
plot(V,eff);
legend('b = 8000 N/(m/s)^.5','b = 10000 N/(m/s)^.5','b = 12000 N/(m/s)^.5','b = 14000 N/(m/s)^.5','b = 16000 N/(m/s)^.5');
ylabel('Efficiency');

title(['Efficiences for Rectified Case at Various Damping Relationships - Dram = ' num2str(D*39.4,2) 'in  Drod = ' num2str(Dr*39.4,2) 'in   PipeSize = ' num2str(pipe_size) 'in']);

%Linear cases
clear eff PressAtMotor PowAtMotor;

PressIdeal = b*V/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);

eff = PowAtMotor./PowIdeal;
ah2 = subplot(3,1,2)
plot(V,eff);
legend('b = 8000 N/(m/s)','b = 10000 N/(m/s)','b = 12000 N/(m/s)','b = 14000 N/(m/s)','b = 16000 N/(m/s)');
ylabel('Efficiency');



%Specific Motor RexRoth-M33
b = 10000;
PressIdeal = b*V/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);
eff = PowAtMotor./PowIdeal;

PP = PressAtMotor/6894;  %Psi
QQ = V*A*(39.4^3);

for i = 1:length(eff)
  [N(i),T(i)] = RexRoth_MR33(PP(i),QQ(i));
end

eff_m = (T.*N)./(6894.757*(.0254^3)*PP.*QQ);

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
plot(V,eff_m,V,eff,V,eff_m.*eff);
xlabel('Ram Speed (m/s)');
ylabel('Efficiency');
title('Force-Speed to Shaft power efficiency for Rexroth-M33 Bore diam = 2.5 in, ram diam = 2 in');



%Specific Motor Parker-TB36
b = 10000;
PressIdeal = b*V/A;
PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);
eff = PowAtMotor./PowIdeal;

PP = PressAtMotor/6894;  %Psi
QQ = V*A*(39.4^3);

for i = 1:length(eff)
  [N(i),T(i)] = Parker_TB36(PP(i),QQ(i));
end

eff_m = (T.*N)./(6894.757*(.0254^3)*PP.*QQ);

figure('Position',[520   142   843   956],'PaperOrientation','portrait','PaperPositionMode','auto');
plot(V,eff_m,V,eff,V,eff_m.*eff);
xlabel('Ram Speed (m/s)');
ylabel('Efficiency');
title('Force-Speed to Shaft power efficiency for Parker-TB36 Bore diam = 2.5 in, ram diam = 2 in');




%Quadratic cases
%clear eff PressAtMotor PowAtMotor PressIdeal PowIdeal F;

%F = b*V.^2;

%PressIdeal = F/A;
%PressAtMotor = PressIdeal-ones(length(b),1)*deltaP;
%PowIdeal = PressIdeal.*(ones(length(b),1)*Q);
%PowAtMotor = PressAtMotor.*(ones(length(b),1)*Q);
%eff = PowAtMotor./PowIdeal;
%ah3 = subplot(3,1,3)
%plot(V,eff);
%legend('b = 8000 N/(m/s)^2','b = 10000 N/(m/s)^2','b = 12000 N/(m/s)^2','b = 14000 N/(m/s)^2','b = 16000 N/(m/s)^2');
%ylabel('Efficiency');
%xlabel('Ram Speed (m/s)');

%linkaxes([ah1 ah2 ah3],'x');