From f7458b7ff764b9606873c764f405a0842ff07683 Mon Sep 17 00:00:00 2001
From: Ildeberto de los Santos Ruiz
 <44104319+isantosruiz@users.noreply.github.com>
Date: Thu, 25 Jul 2019 11:52:32 -0500
Subject: [PATCH] Add files via upload

---
 roughness.m | 85 +++++++++++++++++++++++++++++++++++------------------
 1 file changed, 56 insertions(+), 29 deletions(-)

diff --git a/roughness.m b/roughness.m
index e555094..0b17922 100644
--- a/roughness.m
+++ b/roughness.m
@@ -1,36 +1,63 @@
-function epsilon = roughness(Re,f)
-%ROUGHNESS  Compute the relative roughness coefficient
-%   of a pipe from values of the friction factor and
-%   Reynolds number for different flow conditions.
-%---------------------------------------------------------
-%   Sintax
+function epsilon = roughness(varargin)
+%ROUGHNESS  Compute the relative roughness coefficient of a pipe from
+%   values of the friction factor and Reynolds number for different
+%   operating points.
+%-----------------------------------------------------------
+%   Syntax:
 %      epsilon = roughness(Re,f)
-%---------------------------------------------------------
-%   Arguments
-%           Re : Array of Reynolds numbers.
-%            f : Array of Darcy-Weisbach friction factors.
-%      epsilon : Relative roughness coefficient, k/D.
-%---------------------------------------------------------
-%   Examples
+%      epsilon = roughness(Hin,Hout,Q,nu,L,D,g)
+%-----------------------------------------------------------
+%           Re : Reynolds number
+%            f : Darcy-Weisbach friction factor
+%          Hin : Piezometric head at pipe inlet [m]
+%         Hout : Piezometric head at pipe outlet [m]
+%            Q : Flow rate [m^3/s]
+%           nu : Kinematic viscosity [m^2/s]
+%            L : Pipe length [m] - scalar value
+%            D : Pipe diameter [m] - scalar value
+%            g : Gravity acceleration [m/s^2] - scalar value
+%      epsilon : Relative roughness coefficient
+%-----------------------------------------------------------
+%   Example:
 %      Re = [47525, 74725, 99490, 123013];
 %      f = [0.022786, 0.021086, 0.020241, 0.019698]; 
 %      epsilon = roughness(Re,f)
-%---------------------------------------------------------
-%   (c) 2018, Ildeberto de los Santos Ruiz
-%---------------------------------------------------------
+%-----------------------------------------------------------
+%   Author:
+%      Ildeberto de los Santos Ruiz
+%      idelossantos@ittg.edu.mx
+%-----------------------------------------------------------
 %   Cite as:
-%
-%   Ildeberto de los Santos Ruiz. (2018, November 9).
-%   Friction and Roughness v2.0 (Version v2.0). Zenodo.
-%   http://doi.org/10.5281/zenodo.1481992
-%---------------------------------------------------------
+%      Santos-Ruiz, Ildeberto. (2018, November 9).
+%      Friction and Roughness. Zenodo.
+%      http://doi.org/10.5281/zenodo.1481992
+%-----------------------------------------------------------
+
+switch nargin
+    case 2
+        % roughness(Re,f)
+        Re = varargin{1};
+        f = varargin{2};
+    case 7
+        % roughness(Hin,Hout,Q,nu,L,D,g)
+        Hin = varargin{1};
+        Hout = varargin{2};
+        Q = varargin{3};
+        nu = varargin{4};
+        L = varargin{5};
+        D = varargin{6};
+        g = varargin{7};
+        Re = 4*Q./(pi*nu*D);
+        f = g*pi^2*D^5*(Hin-Hout)./(8*L*Q.^2);
+    otherwise
+        error('Only 2 or 7 input arguments are accepted!')
+        return
+end
+
+colebrook = @(f,Re,epsilon) 1./sqrt(f)+2*log10(epsilon/3.7+...
+    2.51./(Re.*sqrt(f)));
+objfunc = @(epsilon) colebrook(f(:),Re(:),epsilon);
+options = optimoptions('lsqnonlin','Display','none');
+epsilon = lsqnonlin(objfunc,eps,0,1,options);
 
-objfunc = @(epsilon)(LambertW(Re(:),epsilon)-f(:))'*(LambertW(Re(:),epsilon)-f(:));
-options = optimoptions('fsolve','Algorithm','Levenberg-Marquardt','Display','off');
-epsilon = fsolve(objfunc,eps,options);
-    function f = LambertW(Re,epsilon)
-        a = 2.51./Re;
-        b = epsilon/3.7;
-        f = 1./(2*lambertw(0,log(10)./(2*a).*10.^(b./(2*a)))/log(10)-b./a).^2;
-    end
 end
\ No newline at end of file