/*
** program is based on finding roots of the semi-empirical Colebrook
** equation for flow in a circular pipe. this equation is solved iteratively
** using Newton/Raphson method for the diameter pipe required.
**
** user inputs are by changing the #define values of
** flow rate Q required
** length of fuel line L
** tolerable pressure drop over fuel line of delta_p
**
** NOTE: INPUT VALUES AS DECIMAL ie (3.0) NOT (3), AS PROGRAM SOMETIMES NOT HAPPY
** WITH THIS
*/

#include <stdio.h>
#include <math.h>

/* USER INPUTS VIA #DEFINE VALUES BELOW 
*/

/* Q is required flow rate L/min
*/

#define Q (3.1) /* L/min */

/* L is length of fuel line in m
*/

#define L (5.0)	/* m */

/* delta_p is the allowable pressure drop over the fuel line length,
** in psi. ie if p1 is pump output and p2 is pressure
** required at end of fuel line, delta_p=p1-p2
*/

#define delta_p (5.0) /* ie no more than 5psi pressure drop is tolerable */

/* SHOULDN'T HAVE TO CHANGE ANYTHING ELSE BELOW HERE
*/

/* rho is density of fuel kg/m^3
*/

#define rho (680) /* kg/m^3 */

/* mu is absolute viscosity of fuel Ns/m^2
*/

#define mu (2.92e-4)

/* e is average roughness of commercial pipe
e = .046 mm for commercial steel, e = .0015 mm for drawn tubing
*/

#define e (.0015/1000) 	/* metres */


#define PI (3.1415926)

#define Q1 (Q/1000/60) /* L/min converted to m^3/s */
#define del_p (delta_p*6.895e3) /* psi converted to Pa */
#define A (del_p*pow(PI,2)/(8*L*rho*pow(Q1,2)))
#define B (rho*4*Q1/mu/PI)
#define C (pow(A, -.5))
#define D (e/3.7)
#define E (2.51/B*pow(A, -.5))
#define ERROR .000000000001

double fn_value(double x);
double fn_deriv(double x);
double calc_fval(void);

void main()
{
	double blah;

	printf("\n***************************************\n");
   printf("The calculated diameter is %.1fmm\n", calc_fval()*1000);
   printf("***************************************\n\n");

   printf("Press Enter to terminate\n");

   if (scanf("%c", &blah)=='\n');
}

/* fn_value(): returns value of the Colebrook function in terms of pipe
** diameter; with reynolds number, friction factor velocity etc etc
** eliminated
*/

double fn_value(double x)
{
   double fn;

  	fn=C*pow(x,-2.5) + 2*log10(D*pow(x,-1) + E*pow(x,-1.5));

   return fn;
}

/* fn_deriv(): returns the value of the derivative of the Colebrook eqn.
*/

double fn_deriv(double x)
{
   double deriv;

   deriv = -2.5*C*pow(x,-3.5)+2/log(10)*1/(D*pow(x,-1)+E*pow(x,-1.5))*
   		  (-D*pow(x,-2)-3/2*E*pow(x,-2.5));

   return deriv;
}

/* calc_fval calculates and returns the diameter required iteratively
** using the Newton/Raphson scheme.
*/

double calc_fval(void)
{
   double fi, fii;  /* fi = first guess, fii = next guess */

   fi = 1e-4;           /* initial guess of .1mm diameter */

   while (1) {             /* Newton/Raphson scheme */

      fii = fi - fn_value(fi)/fn_deriv(fi);

      if (fabs(fii-fi)/fabs(fii) < ERROR)   /* iteration stop criteria */
         break;
      else
         fi=fii;
   }
   return fii;
}