"""
Demostration of the Hardy-Cross relaxation solution of the pipeline network
(1) (2) (3) (7)
----o-----o-----o------
\(4) |(5) /(6)
\ | /
\--o--/
"""
import numpy as n
import pylab as p
L = n.array([2,2,3,3,2,3,2],dtype='f')
D = n.array([25,25,20,20,15,20,25],dtype='f')
C = n.array([120,100,110,120,130,120,130],dtype='f')
r = lambda l,d,c: 1.526e7/(c**1.852*d**4.87)*l
hf = lambda R,Q: R*n.sign(Q)*n.abs(Q)**1.852
R = r(L,D,C)
branch = n.array([[2,5,4],[3,6,5]])-1
rows,cols = branch.shape
Q = n.array([500,250,250,-250,0,-250,500],dtype='f')
dQ = 500.0
while abs(dQ) > 0.5:
for i in n.arange(rows):
y = hf(R,Q)
yq = n.abs(1.852*y/n.abs(Q))
yq[n.isnan(yq)] = 0.0
sumyq = n.sum(yq[branch[i,:]])
sumy = n.sum(y[branch[i,:]])
dQ = -1*sumy/sumyq
print("dQ = %f" % dQ)
Q[branch[i,:]] += dQ
print("Discharges [m^3/hr]:\n")
print Q
Deq = 25
Ceq = 120
y = hf(R[1],Q[1])+hf(R[2],Q[2])
Leq = y/(r(1,25,120)*Q[0]**1.852)
Leq = Leq + L[0] + L[6]*(120/C[6])**1.852
print("Equivalent length = %3.2f km" % Leq)