Programando en Octave y C - Diagrama de Moody

Estimados, Después de publicar temas relacionados a electrónica, informática, redes, etc; hoy presento el desarrollo de un tema totalmente diferente. La aplicación de la programación en C y Octave para solucionar problemas de Ingeniería Hidráulica. Esto a raíz de que mi padre está cursando su Maestría en Ingeniería Hidráulica en la UNI - Perú y necesitaba una herramienta que le permita calcular lo que a continuación expongo. En el curso de flujos en superficies libres tocaron el tema del uso del Diagrama de Moody para calcular el factor de fricción de un determinado fluido bajo ciertas condiciones. Esto podía tener dos escenarios dependiendo de los valores obtenidos preliminarmente. Podía ser un escenario de flujo laminar o flujo turbulento y ésto depende del rango en que se encuentre el número de Reynolds. Si éste número se encuentra entre 0 y 2300 e considera que el fluido se encuentra en régimen laminar y la fórmula para obtener el factor de fricción es : f = 64/nReynolds. Si está por encima de los 2300, se considera que el fluido se encuentra en régimen turbulento y se utiliza la fórmula de Colebrook para obtener el factor de fricción, además de tener en cuenta otra variable conocida como Rugosidad relativa (ed):

La cual es necesario iterar hasta encontrar el factor.
Alternativamente se puede usar la expresión propuesta por Haaland, en el que logra despejar "f" de la ecuación de Colebrook.

He usado Octave (altamente compatible con matlab) para generar el diagrama de Moody.
Lo interesante con esto, es que aparte de ayudar a mi padre, he aprendido comandos interesantes como num2str, loglog, text y he reforzado mi conocimiento de bucles anidados y uso de arrays en Octave/Matlab.
Expongo el código por si alguien esté interesado.
Con ese código se puede obtener el diagrama de Moody para régimen turbulento, el cual se puede usar para hallar aproximadamente el factor de fricción para un determinado número de Reynolds y rugosidad relativa.

%%PROGRAMA QUE GRAFICA EL DIAGRAMA DE MOODY - FLUJO TURBULENTO
%%SE ITERA CON DETERMINADOS VALORES DE RUGOSIDAD RELATIVA E/D vs nReynolds
%%EL NUMERO DE REYNOLDS SE CALCULA EN BASE A LA FORMULA:
%%     nReynolds = VELOCIDAD*DENSIDAD*DIAMETRO_TUBERIA/MU
%% LA RUGOSIDAD RELATIVA SE CALCULA EN BASE A LA FORMULA:
%%     e/DIAMETRO = ed
%% EL COEFICIENTE DE RUGOSIDAD (f) SE CALCULA EN BASE A LA FORMULA DE HAALAND:
%% 1/f^(1/2)=-1.8*log[(6.9/nReynolds)+(ed/3.7)^1.11]
clear all;close all;clc;
nReynolds = [3000];
ed = [0.00001 0.00002 0.00005 0.0001 0.0002 0.0004 0.0006 0.0008 0.001 ...
0.0015 0.002 0.003 0.004 0.006 0.008 0.01 0.0125 0.015 0.0175 0.02 0.025 ...
0.03 0.035 0.04 0.045 0.05 0.06 0.07];
f = [];
for i = 2:214
  nReynolds = [nReynolds 1.05015*nReynolds(i-1)];
end
hold on;
for i = 1:length(ed)
  for j = 1:length(nReynolds)
  f = 0.30864/(log10((6.9/nReynolds(j))+(ed(i)/3.7)^1.11))^2;
  eval(num2str(i,'f%d(j)=f;'));
  end
end
figure (1)
hold on;grid on;
for i=1:length(ed)
eval(num2str(i,'loglog(nReynolds,f%d)'))
end
title({'Diagrama de Moody - Seccion Flujo Turbulento'; ' '; ...
'{f} ^{-1/2} = -1.8 * log10 [ ( 6.9 / nReynolds ) + ( ed / 3.7 ) ^{1.11}]'})
xlabel ('Number Reynolds'); ylabel ('Factor friction'); 
for i=1:length(ed)
VAL = ed(i);
str = ['ed = ',num2str(VAL)];
text(nReynolds(214),eval(num2str(i,'f%d(214)')),str);
end

Con este algoritmo, se obtiene la siguiente gráfica:

Si necesitamos un programa que calcule analíticamente el factor de fricción para ambos regímenes, presento este programa hecho en C y compilado con gcc bajo Debian 8.2.

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define DENSIDAD 1.23
#define MHU 0.0000179

void muestraEncabezado ()
{
puts("UNIVERSIDAD NACIONAL DE INGENIERIA\t\t MAESTRIA EN INGENIERIA HIDRAULICA\n");
    puts("\n\t\t\t PROGRAMA QUE DETERMINA EL TIPO \n\t\t  DE REGIMEN DEL AIRE A TRAVES DE UNA TUBERIA \
    \n\t\t\tY CALCULA EL FACTOR DE FRICCIÓN\n");
    puts("Alumno    : Jorge Piscoya Fernández");
    puts("Curso     : Flujo en Superficie Libre");
    puts("-------------------------------------\n");
}
float mostrarReynolds(float diametro, float velocidad)
{
    return (diametro*velocidad*DENSIDAD)/(MHU);
}
float regimenLaminar (float diametro, float velocidad)
{
    return 64 / mostrarReynolds(diametro, velocidad);
}
float regimenTurbulento (float ed, float diametro, float velocidad)
{
    float f;
    float nR;
    nR = mostrarReynolds(diametro, velocidad);
    f = 0.30864/(pow(log10((6.9/nR)+(pow(ed/3.7,1.11))),2));
    return f;
}
int main()
{
    float diametro = 0;
    float velocidad = 0;
    float nReynolds = 0;
    muestraEncabezado();
    printf("Ingrese diámetro de tubería en (m): ");
    scanf("%f", &diametro);
    printf("\n");
    printf("Ingrese la velocidad del fluido (m/s): ");
    scanf("%f", &velocidad);
    printf("\n");
    nReynolds = mostrarReynolds(diametro, velocidad);
    printf("El número de Reynolds para el fluido ingresado es: %f\n", nReynolds);
    if (nReynolds < 2300)
    {
        puts("El fluido se encuentra en régimen laminar");
        printf("El coeficiente de rugosidad es: %f\n", regimenLaminar(diametro, velocidad));
    }
    else
    {
        float ed, f;
        puts("El fluido se encuentra en régimen turbulento\n\n");
        puts("Se procederá a calcular el factor de fricción en Reg. Turbulento\n");
        printf("Ingrese el valor de la Rugosidad Relativa (ed) del medio: ");
        scanf("%f", &ed);
        f = regimenTurbulento(ed, diametro, velocidad);
        printf("\nEl factor de fricción en regimen turbulento es: %f\n",f);
    }
    puts("\nPresione ENTER para salir\n");
    while(getchar() != '\n');
    getchar();
    return 0;
}

Y se obtiene la siguiente ventana:

Voy a intentar meter ncurses para que sea más colorida la interface, pero será en otro momento.
Para ejecutarlo en Windows, pueden usar DEV-C++ ya que es parte del standar ANSI C, así podrán obtener el ejecutable para windows.

Para poder realizar este trabajo, he seguido este manual teórico más la explicación de mi padre:
http://es.scribd.com/doc/67350243/Diagrama-de-Moody

Con esto reforzado mi conocimiento básico de manejo de funciones definidas para el desarrollo del algoritmo. Si tienen alguna sugerencia de cómo mejorar el código, es bienvenido. Saludos