"""
Implementation of the Hardy-Cross relaxation solution of the pipeline network,
given in the following example:
(2)
A------------B
| / |
| / |
| / |
|(1) (3)/ (4)|
| / |
| / |
| / |
| / |
| / |
| / (5) |
C------------D
QA = 500 # m**3/h
QB = 0 # m**3/h
QC = 0 # m**3/h
QD = -500 # m**3/h
R1 = 0.000298
R2 = 0.000939
R3 = 0.00695
R4 = 0.000349
R5 = 0.000298
"""
from numpy import *
QA = 500
QB = 0
QC = 0
QD = -500
L = array([2,3,3,2,3],dtype='f')
D = array([25,20,20,15,20],dtype='f')
C = array([100,110,120,130,120],dtype='f')
r = lambda L,D,C: L*1.526e7/(C**1.852*D**4.87)
R = r(L,D,C)
print "Resistances: ", R
hf = lambda R,Q: R*sign(Q)*power(abs(Q),1.852)
branch = array([[2,3,1],[3,4,5]])-1
rows,cols = branch.shape
Q = array([-250, 250.0, 0.0, 250, -250.0])
dQ = 1.0
while abs(dQ) > 0.5:
for i in arange(rows):
y = hf(R,Q)
yq = abs(1.852*y/abs(Q))
yq[isnan(yq)] = 0.0
sumyq = sum(yq[branch[i,:]])
sumy = sum(y[branch[i,:]])
dQ = -1*sumy/sumyq
print("dQ = %f" % dQ)
Q[branch[i,:]] += dQ
print("Discharges [m^3/hr]:\n")
print Q