function Prob_5_2c
clear, clc, format short g, format compact
D_list=[4.026/12 5.047/12 6.065/12 7.981/12]; % Inside diameter of pipe (ft)
T = 60; %Temperature (deg. F) 
for i = 1:4
    D = D_list(i);
    j=0;
    for L=500:500:10000
        j = j+1;
        L_list(j)=L; % Effective length of pipe (ft)
        [v(j,i),fval]=fzero(@NLEfun,[1 20],[],D,L,T);
        if abs(fval)>1e-10
            disp([' No Convergence for L = ' num2str(L) ' and D = ' num2str(D)]);
        end
        q(j,i) = v(j,i) * pi * D ^ 2 / 4* 7.481 * 60; %Flow rate (gpm) 
    end
end
disp('                 Flow Velocity (ft/s) versus Pipe Length and Diameter');
disp('                                   Tabular Results');
disp('');
disp('         L          D=4"         D=5"         D=6"          D=8"');
Res=[L_list' v];
disp(Res);
plot(L_list,v(:,1),'-',L_list,v(:,2),'+',L_list,v(:,3),'*',L_list,v(:,4),'x');
legend(' D=4"',' D=5"',' D=6"',' D=8"');
title(' Flow Velocity')
xlabel('Pipe Length (ft)');
ylabel('Velocity (ft/s)');
pause
disp('                    Flow Rate (gpm) versus Pipe Length and Diameter');
disp('                                   Tabular Results');
disp('');
disp('         L\D        D=4"         D=5"         D=6"         D=8"');
Res=[L_list' q(:,1) q(:,2) q(:,3) q(:,4)];
disp(Res);
plot(L_list,q(:,1),'-',L_list,q(:,2),'+',L_list,q(:,3),'*',L_list,q(:,4),'x');
legend(' D=4"',' D=5"',' D=6"',' D=8"');
title(' Flow rate')
xlabel('Pipe Length (ft)');
ylabel('Flow rate (gpm)');
function fv = NLEfun(v,D,L,T)
epsilon = 0.00015;%Surface rougness of the pipe (ft) 
rho = 62.122 + T * (0.0122 + T * (-0.000154 + T * (0.000000265 - (T * 0.000000000224)))); %Fluid density (lb/cu. ft.) 
deltaz = 300; %Elevation difference (ft) 
deltaP = -150; %Pressure difference (psi)
vis = exp(-11.0318 + 1057.51 / (T + 214.624)); %Fluid viscosity (lbm/ft-s) 
pi = 3.1416; %The constant pi 
eoD = epsilon / D; %Pipe roughness to diameter ratio (dimensionless) 
Re = D * v * rho / vis; %Reynolds number (dimesionless) 
if (Re < 2100) %Fanning friction factor (dimensionless) 
    fF = 16 / Re;
else
    fF = 1 / (16 * log10(eoD / 3.7 - (5.02 * log10(eoD / 3.7 + 14.5 / Re) / Re)) ^ 2);
end 
fv = v - sqrt((32.174 * deltaz + deltaP * 144 * 32.174 / rho) / (0.5 - (2 * fF * L / D))); %Flow velocity (ft/s)